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SUMMARY 


The  objective  of  the  research  in  this  program  was  to  develop  Langrangian  meth¬ 
ods  on  triangular  grids  and  apply  these  methods  to  the  calculation  of  life-history  and 
dynamics  of  fuel  droplets.  During  the  research  period,  major  advances  were  made  both 
in  numerical  technology  and  in  the  solution  of  problems  so  difficult  that  they  have  not 
been  attempted  before. 

With  respect  to  numerical  technology,  the  two-dimensional  code  SPLISH  was  con¬ 
verted  to  a  VAX  and  then  to  a  CRAY  computer.  New  graphics  systems  were  developed. 
Further  testing  of  the  basic  SPLISH  hydrodynamic  algorithms  as  well  as  the  surface 
tension  algorithm  were  performed  on  internal  gravity  and  capillary  waves.  A  reor¬ 
ganization  of  the  computer  code  itself  is  nearly  complete.  This  will  make  the  code 
user-friendly  and  portable.  Now  it  should  be  much  easier  to  use,  and  therefore  useable 
on  many  new  kinds  of  problems. 

First,  previously  calculated  flows  of  the  distortion  and  breakup  of  a  droplet  due 
to  differences  in  flow  velocities  between  the  droplet  and  the  external  media  were  re¬ 
computed  to  verify  the  conversion.  Then  a  number  of  calculations  of  droplet  distortion 
and  breakup  due  to  shear  flows  were  made.  Qualitative  comparisons  to  experimental 
results  were  made  for  the  case  when  the  droplet  density  and  external  fluid  density  were 
nearly  equal.  Our  calculation  and  the  experiments  by  Mason  and  coworkers  showed 
small  droplets  torn  off  the  large  drop  by  the  forces  in  the  shear  flow.  A  preliminary 
calculation  of  a  droplet-droplet  collision  shows  the  distortion  of  droplets  before  they 
collide.  Forced-flow  and  inflow-outflow  boundary  conditions,  needed  to  do  quantitative 
comparisons  to  experimental  shear  flows,  were  added  to  the  model. 

RESEARCH  OBJECTIVES 

Our  major  research  objective  was  to  develop  models  that  could  be  used  to  give 
the  basic  information  needed  for  constructing  appropriate  models  for  dense  sprays. 
This  is  necessary  because  the  limiting  assumptions  in  current  drop  models  imply  that 


the  predictions  given  by  spray  models  in  dense  or  high-pressure  spray  regimes  are 
not  as  accurate  as  they  could  be.  Better  drop  models  should  be  formulated  using 
information  about  collision  and  breakup  rates  relating  to  droplet  size  and  velocity. 
However  such  rate  calculations  are  complicated  by  a  rich  variety  of  breakup,  collision, 
and  coalescence  modes  that  can  occur  in  different  pressure  and  velocity  regimes.  For 
example,  collisions  can  result  in  coalescence,  elastic  rebound,  inelastic  scattering  with 
increased  drop  deformation  and  oscillation  or  shattering  into  two  or  more  droplets. 
In  addition,  rapidly  rotating  or  highly  deformed  drops  may  split  into  several  drops, 
shattering  may  occur  at  critical  pressures,  and  hydrodynamic  breakup  of  drops  may 
occur  through  several  modes,  including  the  bag  breakup  mode  and  droplet  stripping 
due  to  velocity  gradients  in  the  exterior  flow  field.  Drop  transport  is  directly  affected 
by  these  phenomena,  and  the  presence  of  the  drops  affects  the  flow  field  in  both  the 
mean  flow  and  turbulent  fluctuations. 

Our  approach  was  to  develop  the  numerical  technology  associated  with  Lagrangian 
methods  on  triangular  grids  to  the  point  where  it  could  be  used  to  model  the  life-history 
and  dynamics  of  fuel  droplets.  In  particular  these  methods  would  be  used  to  determine 
the  conditions  under  which  fuel  droplets  greatly  distort  and  shatter,  as  would  occur  in 
shear  flows  and  droplet  collisions. 

WORK  PERFORMED 

The  basis  of  the  numerical  method  is  a  Lagrangian  convection  algorithm  which 
uses  a  triangular  grid  instead  of  a  quadrilateral  grid.  In  this  method,  the  grid  dy¬ 
namically  restructures  itself  according  to  preprogrammed  criteria.  In  addition,  new 
triangles  are  added  or  old  triangles  are  deleted  to  change  the  resolution,  also  according 
to  preprogrammed  criteria.  Triangle  deletion  allows  the  evolution  of  singly  connected 
regions  to  multiconnected  regions,  such  as  would  occur  when  droplets  break  up,  or  the 
evolution  of  multiconnected  regions  to  singly  connected  regions,  such  as  would  occur 
during  droplet  collision.  Since  the  algorithm  is  Lagrangian,  interfaces  can  be  tracked 
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expicitly.  Surface  tension  forces  and  viscous  forces  are  included  in  addition  to  the 
Lagrangian  convection. 

Conversion  of  the  Triangular  Langrangian  Code  SPLISH  to  Other  Computers 

In  order  to  begin  the  work  on  this  project,  we  had  to  move  the  code  SPLISH 
originally,  written  for  the  TI-ASC,  first  to  a  VAX  and  then  to  the  CRAY  computer. 
In  addition  to  converting  the  code,  we  had  to  develop  different  graphics  packages  for 
the  new  computers.  This  was  an  extensive  effort.  Along  the  way,  many  of  the  original 
calculations  done  under  NASA  support  were  redone  and  checked.  In  particular,  the 
droplet  oscillation  test  calculations  and  one  of  the  calculations  of  a  flow  pjist  a  droplet 
were  performed.  These  calculations  are  reported  in  detail  in  Appendix  B,  which  is  a 
copy  of  the  paper  to  be  published  in  1988  in  the  Journal  of  Computational  Physics, 
and  to  a  lesser  extent  in  Appendix  A,  which  is  paper  number  AI.\A-87-0539  from  the 
AIAA  Aerospace  Sciences  meeting  of  January,  1987. 

Detailed  Quantitative  Tests  of  the  Surface  Tension  Algorithm 

Developing  an  accurrate  algorithm  to  model  the  effects  of  interfacial  surface  tension 
forces  was  much  more  difficult  than  we  expected.  We  believe  that  this  had  never 
been  satisfactorily  addressed  in  the  literature.  We  tested  many  different  approaches, 
including  the  one  we  finally  used,  as  documented  in  Appendix  B.  At  first  we  used  an 
n  =  2  normal  mode  droplet  oscillation  as  a  test  of  the  accuracy  of  the  surface  tension 
algorithm.  The  algorithm  we  finally  used  was  accurate  enough  to  allow  us  to  compute 
a  number  of  periods  of  oscillation  of  the  droplet.  However,  the  convergence  properties 
of  the  algorithm  were  tested  by  comparing  theoretical  and  computational  calculations 
of  the  properties  of  internal  capillary  waves.  The  result  was  that  the  surface  tension 
algorithm  is  first  order  in  the  mesh  size,  although  the  convection  algorithm  is  second 
order. 
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There  is  one  other  result  implied  by  our  tests  of  the  surface  tension  algorithm. 
The  algorithm  works  well  when  the  surfaces  are  not  too  highly  distorted.  For  highly 
distorted  surfaces  the  time  and  length  scales  associated  with  the  surface  tension  phe¬ 
nomenon  are  much  smaller  than  the  convective  time  and  length  scales.  For  droplet 
breakup  the  physics  of  interest  happens  on  the  convective  scales.  For  distorted  surfaces 
many  more  points  and  time  steps  are  needed  to  compute  the  motion  due  to  the  surface 
tension  forces  accurately,  than  are  needed  to  compute  the  convective  motion  accurately. 

The  details  of  this  work  can  be  found  in  Appendix  A.  These  results  have  been 
presented  at  the  American  Physical  Society  Meeting  of  the  Division  of  Fluid  Dynam¬ 
ics,  November,  1985,  and  at  the  International  Symposium  on  Computational  Fluid 
Dynamics,  Tokyo,  in  August,  1985. 

Calculations  of  Droplet  Distortion  and  Breakup  in  Shear  Flows 

Calculations  of  silicon  oil  droplets  were  compared  to  the  analytic  and  experimental 
work  of  Mason  and  co-workers.  In  these  calculations,  a  drop  of  density  0.98  g/cc  and 
diameter  1  mm  was  placed  in  an  initial  shear  flow  prescribed  by  =  G(y  -  yj),  for 
points  outside  the  droplet  and  =  0  for  points  inside  the  droplet.  The  parameter  yj 
is  is  the  y-coordinate  of  the  center  of  the  droplet  and  G  gives  the  magnitude  of  the 
shear.  The  y-coordinate  of  the  velocity  is  initially  zero  everywhere. 

Figure  1  shows  the  time  evolution  of  one  such  calculation  in  which  the  viscosity 
and  small  surface  tension  coefficients  were  small.  The  experimental  results  for  this  case 
showed  that  the  droplet  stretches  in  the  direction  of  the  flow  and  small  droplets  break 
from  the  tips.  This  is  similar  to  the  calculated  results  shown  in  Figure  1.  The  small 
droplets  shed  in  the  calculation  were  larger  than  those  shed  in  the  experiments.  The 
discrepancy  is  due  to  the  resolution  of  the  numerical  calculation,  in  which  we  have  only 
a  few  triangles  available  to  represent  the  shed  droplets.  More  resolved  calculations, 
which  we  are  currently  doing,  should  provide  a  more  solid  basis  for  comparison  with 
the  experimants. 
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FIGURE  1 


New  and  unique  calculations  of  the  breakup  of  a  dense  fuel  droplet  in  a  hot  air  shear 
flow  were  also  performed.  Here  the  stretching  of  the  drop  in  the  direction  of  the  shear 
is  much  less  pronounced  due  to  the  large  density  of  the  drop  relative  to  the  background 
air.  Better  resolution  has  allowed  us  to  see  smaller  droplets  being  shed.  Currently,  very 
resolved  calculations  are  being  performed.  The  results  of  these  calculations  have  been 
reported  at  the  1986  SIAM  National  Meeting,  the  American  Physical  Society  Meeting  of 
the  Division  of  Fluid  Dynamics,  November,  1986,  and  the  Aerospace  Sciences  meeting 
of  the  American  Institute  of  Aeronautics  and  Astronautics,  January,  1987.  A  detailed 
discussion  of  one  of  these  calculations  can  be  found  in  Appendix  A. 

Calculations  of  Droplet -Droplet  Collisions 

A  preliminary  calculation  of  a  droplet-droplet  collision  was  performed.  This  prob¬ 
lem  consists  of  a  stationary  target  droplet  and  a  moving  projectile  droplet.  In  the 
calculation,  the  two  kerosene  droplets  were  placed  in  a  background  fluid  of  hot  air. 
The  target  droplet  is  initially  at  rest,  while  the  background  fluid  and  the  projectile 
droplet  on  the  left  were  initialized  with  a  potential  flow  of  1  m/s  about  the  target 
droplet.  Figure  2  shows  the  initial  phases  of  this  calculation.  As  the  projectile  nears 
the  target,  both  droplets  deform  and  present  nearly  flat  faces  to  each  other.  The  lower 
parts  of  the  droplets  begin  to  merge  first.  This  merger  is  accomplished  numerically  by 
the  deletion  of  small  triangles  of  background  fluid  between  the  two  droplet  interfaces. 
Since  the  new  interface  now  has  a  drastic  change  in  curvature,  the  surface  tension  forces 
try  to  expel  the  background  fluid  upwardly  from  between  the  two  main  sections  of  the 
kerosene  mass.  This  process  produced  large  numerically  generated  fluctuations  in  the 
pressures  in  that  region. 

Analyzing  why  this  happened  provides  important  information  about  problems  in 
the  basic  algorithm  and  potential  solutions  of  these  problems.  First,  there  are  potential 
errors  in  the  procedure  used  to  delete  triangles  near  interfaces.  These  would  result  in 
a  phase  lag  in  the  pressures  at  newly  created  interfaces.  A  second  possibility  has  its 


origin  in  the  values  of  the  triangle  velocities  after  vertex  deletion.  The  triangle  velocities 
determine  the  rate  of  expansion  or  contraction  of  the  volumes  about  each  vertex.  The 
pressures  react  to  prevent  this  change.  When  a  triangle  is  deleted,  the  velocities  of 
the  neighboring  triangles  must  be  altered  to  maintain  conservation  of  momentum  and 
circulation.  If  these  velocities  are  changed  in  such  a  way  as  to  produce  a  huge  volume 
expansion,  we  would  obtain  the  symptoms  we  observe. 

These  possibilites  are  being  evaluated,  several  proposed  fixes  are  being  tested,  and 
the  calculation  will  proceed  once  we  have  an  answer  to  the  problem. 

Addition  of  Boundary  Conditions 

The  experiments  of  Mason  and  coworkers  of  droplet  shears  were  performed  pri¬ 
marily  with  highly  viscous  silicon  oils  with  the  shear  flow  maintained  by  the  motion 
of  the  containment  device.  Since  the  visous  terms  are  sinks  for  energy  in  the  SPLISH 
algorithm,  highly  viscous  flows  damp  very  quickly  to  no  flow  at  all.  Therefore  boundary 
conditions  which  simulate  the  motion  of  the  experimental  container  have  been  added 
to  the  SPLISH  algorithm  in  order  to  compare  the  code  results  to  the  experiments. 
To  implement  this  boundary  condition  a  layer  of  triangles  is  placed  just  outside  the 
computational  region  on  the  moving  boundary.  These  triangles  have  a  density  equal 
to  that  just  inside  the  moving  boundary  and  are  given  a  prescribed  velocity  parallel 
to  the  boundary.  Through  the  viscous  terms  in  the  momentum  equation,  this  kinetic 
energy  is  then  transferred  to  the  background  fluid.  To  model  the  experiments,  we  now 
set  the  initial  flow  to  zero  everywhere  except  on  the  exterior  boundary  triangles.  The 
shear  flow  experiments  used  impulsive  starts,  which  corresponds  to  giving  the  bound¬ 
ary  layer  of  triangles  a  constant  velocity  initially,  and  ramp  starts,  which  corresponds 
to  giving  the  boundary  layer  of  triangles  a  linear  time  dependent  velocity  up  to  the 
desired  boundary  velocity. 

We  are  also  implementing  inflow-outflow  boundary  conditions.  There  are  problems 
with  doing  this,  some  particular  to  our  formulation,  and  others  because  it  is  simply  a 


very  difficult  problem  (see,  for  example,  Oran  and  Boris,  1987).  One  problem  partic¬ 
ular  to  our  formulation  is  geometric  and  associated  with  the  addition  and  deletion  of 
mesh  points  at  the  boundaries.  Another  is  associated  with  conservation  of  circulation 
and  vertex  cell  area  because  of  the  incomplete  vertex  cells  for  those  vertices  on  the 
boundaries.  A  generic  problem  is  the  physical  boundary  conditions  themselves  which 
determines  how  variables  just  inside  the  boundary,  as  well  as  the  variables  at  newly 
created  vertices,  are  time  advanced.  The  problem  is  made  more  difficult  by  the  irreg¬ 
ular  mesh,  since  there  is  essentially  no  literature  on  inflow-outflow  for  incompressible 
flows  differenced  over  a  general  connectivity  mesh.  Research  on  boundary  conditions 
in  computational  fluid  dynamics  is  an  important  ongoing  area  of  numerical  research. 

The  geometrical  aspects  of  inflow-outflow  boundary  conditions  we  are  now  using 
have  been  tested  and  work  well.  The  algorithm  deletes  those  triangles  with  no  vertex 
within  the  computational  region  and  adds  triangles  at  the  inflow  boundary  whenever 
any  boundary  vertex  is  completely  within  the  computational  region.  The  next  step  is 
both  easier  and  more  difficult:  placing  more  physically  correct  values  of  the  physical 
variables  on  the  new  triangles  and  vertices.  At  the  inflow  boundary,  placing  a  prescribed 
value  of  velocity  on  a  bordering  triangle  may  not  produce  a  vanishing  divergence  on  a 
nearby  interior  vertex.  This  means  that  the  inflow  triangle  velocities  need  to  be  rotated 
slightly  to  produce  the  appropriate  incompressible  flow.  On  the  outflow  boundary,  the 
physical  variables  need  to  be  updated  in  a  manner  that  will  not  produce  numerically 
reflected  waves  off  that  boundary.  The  form  of  the  physical  outflow  boundary  conditions 
that  would  be  best  to  use  in  our  case  are  the  Orlansky  outflow  condition  as  formulated 
by  Chan  in  his  article  Finite  Difference  Simulation  of  the  Planar  Motion  of  a  Ship.  The 
positive  aspect  of  his  approach  is  that  is  upwind  so  there  are  no  reflected  waves  and 
hence  we  do  not  need  to  introduce  artificial  damping. 


Extensions  of  the  SPLISH  Algorithm  to  Three  Dimensions 

Extensions  of  the  Lagrangian  algorithm  on  triangular  grids  to  a  Lagrangian  algo¬ 
rithm  on  tetrahedral  grids  has  always  been  a  major  goal  of  the  work.  However,  as  we 
progressed  into  the  problem,  we  realized  what  a  major  effort  in  algorithm  development 
is  required.  Our  initial  plans  were  that  we  would  have  additional  funding  from  other 
sources  that  we  use  for  algorithm  development,  which,  and  in  the  last  year  of  this  pro¬ 
posal,  would  allow  us  to  calculate  the  properties  of  three-dimensional  droplets  instead 
of  two-dimensional  cylinders. 

The  major  algorithmic  problem  with  Lagrangian  tetrahedra  involve  the  reconnec¬ 
tions.  In  two-dimensions,  the  test  for  reconnection  of  a  side  is  merely  the  comparisons 
of  two  sets  of  opposing  angles  in  the  quadrilateral  formed  by  the  triangles  on  either 
side,  and  the  reconnection  is  merely  the  swap  of  the  diagonals  of  the  quadrilateral.  Re¬ 
connection  in  two  dimensions  is  local  to  the  side  in  question.  In  three  dimensions,  there 
is  no  equivalent  diagonal  swap.  Instead,  one  must  consider  several  tetrahedra  and  faces 
simulateously,  in  order  to  get  the  optimum  grid  connectivity.  Currently,  Dr.  Martin 
Fritts  of  SAIC  has  been  working  on  this  problem  in  collaboration  with  and  with  funding 
from  Lawrence  Livermore  National  Laboratories  and  Los  Alamos  National  Laborato¬ 
ries.  However,  we  at  NRL  are  still  somewhat  skeptical  about  being  able  to  apply  the 
three-dimensional  SPLISH  analog  in  an  efficient  way  in  the  near  future. 

Therefore  we  have  taken  two  other  approaches.  The  first  still  involves  triangles  and 
tetrahedra:  we  have  successfully  developed  a  finite-element  version  of  the  high-order 
monotone  scheme  Flux-Corrected  Transport,  which  uses  triangles  instead  of  quadri¬ 
laterals  as  the  basic  elements.  This  is  bcisically  an  Eulerian  method,  although  a  new 
arbitrary  Lagrangian-Eulerian  (ALE)  now  exists.  This  work  is  being  carried  out  at 
NRL  by  Dr.  Rainald  Lohner  under  the  sponsorship  of  both  D.\RP.A  and  N.\S.\.  These 
methods  now  work  well  and  are  incorporated  in  production  codes.  Combined  with  a 
front-tracking  algorithm,  this  approach  could  now  be  used  to  calcualte  droplet  dynam¬ 
ics  in  three  dimensions. 


The  second  numerical  approach  is  still  in  the  conceptual  stage.  This  involves  using 
the  Monotonic  Lagrangian  Grid  representation,  developed  by  Jay  Boris  at  NRL,  as  a 
basis  for  Lagrangian  fluid  dynamics.  This  is  an  approach  we  hope  to  evaluate  in  the 
next  couple  of  years. 

Continued  Efforts  to  Make  the  Code  More  User-Friendly 

Our  objective  here  is  to  produce  a  computer  code  which  is  more  flexible  in  the 
problems  it  can  handle  as  well  as  more  portable.  The  program  SPLISH  has  the  potential 
to  solve  a  large  class  of  problems.  However,  because  of  the  developmental  nature  of 
the  work,  the  code  became  a  patchwork  of  additions  and  fixes.  In  order  to  apply 
the  algorithm  with  some  facility  to  a  variety  of  problems,  some  efforts  in  structured 
programming  was  required.  The  code  itself  must  be  suitably  structured  so  that  the 
modifications  can  be  made  in  a  straightforward  manner.  .Approximately  80%  of  this 
task  is  complete.  The  only  major  items  remaining  to  be  upgraded  are  the  vertex 
deletion  algorithms  for  small  triangles.  This  will  be  complete  in  FY’88. 

Discussion 

During  the  contract  period,  major  advances  have  been  made  in  this  program  both 
in  numerical  technology  and  in  the  numerical  solution  of  problems  so  difficult  that 
they  have  not  been  attempted  before.  The  physical  problems  are  the  droplet  breakup 
in  shear  flows  and  droplet-droplet  collisions.  The  numerical  problems  involve  accurate 
formulations  of  surface  tension,  boundary  conditions,  and  grid  restructuring  algorithms 
on  the  Lagrangian  triangular  grid.  This  coming  year,  in  FY88,  the  results  of  the 
droplets  in  shear  flows  will  be  written  up  and  submsitted  to  the  Journal  of  Fluid 
Mechanics.  The  calculations  of  droplet-droplet  collisions  will  be  continued  until  the 
numerical  questions  are  answered  and  a  case  is  run. 

While  the  calculations  of  droplet-droplet  collisions  has  only  just  begun,  the 
prospects  for  significant  research  results  is  great.  The  recent  experimental  results  C.K. 
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Law  from  the  University  of  California,  Davis  and  from  N'assar  Ashgriz  from  the  State 
University  of  New  York  at  Buffalo  will  allow  us  some  benchmark  examples  for  further 
code  verification  before  poceeding  to  cases  for  which  there  are  no  experimental  results. 
We  have  already  started  talking  with  Professor  Law  and  we  are  planning  a  case  to 
simulate. 

The  code  reorganization  star^^-d  under  the  sponsorship  of  AFOSR  will  have  the 
most  lasting  benefit.  With  the  capability  of  easy  modification  of  a  computationally 
complex  algorithm,  the  number  of  problems  attacked  by  this  approach  could  multiply. 
One  such  problem  is  the  behavior  of  a  reacting  liquid-gas  interface.  The  computational 
aspects  of  this  problem  is  suited  to  the  SPLISH  algorithm.  Algorithms  for  the  con¬ 
densation  of  the  gas  (or  equivalently  the  evaporation  of  the  liquid)  would  have  to  be 
added.  This  is  a  problem  of  current  interest  to  the  Navy  in  connection  with  liquid-metal 
combustion. 

PERSONNEL 

The  funding  for  this  contract  has  supported  the  research  of  Drs.  David  Fyfe  and 
Elaine  Oran.  Dr.  Martin  Fritts,  who  currently  works  for  Science  Applications,  partic¬ 
ipated  in  the  program  as  a  consultant  during  fiscal  year  1985.  In  addition,  Dr.  J.P. 
Boris  and  Dr.  Rainald  Lohner  regularly  consulted  on  this  work. 

PRESENTATIONS  AND  PUBLICATIONS 

Simulations  of  Two-Dimensional  Fuel  Droplet  Flows.  M.J.  Fritts.  D.E.  Fyfe,  and  E.S. 
Oran.  Prog.  Astro.  Aero.  95.  540-553.  Dynamics  of  Flames  and  Reactive  Systems. 
edited  by  J.R.  Bowen.  N.  .Manson.  .A.K.  Oppenheim.  and  R.l.  Soloukhin.  .ALA.A. 
1984. 

.Application  of  Two-Dimensional  Techniques  to  Combustion,  invited  presentation  at 
Heidelberg  University,  Heidelberg,  West  Germany,  1984. 


Numerical  Simulations  of  Fuel  Droplet  Flows  using  a  Lagrangian  Triangular  Mesh, 
M.J.  Fritts,  D.E.  Oran,  and  E.S.  Oran,  Ninth  International  Conference  on  Numer¬ 
ical  Methods  in  Fluid  Dynamics  pp.  219-224,  Springer- Verlag,  New  York,  1985. 


Multidimensional  Simulation  of  Flames  and  Detonations,  E.S.  Oran,  invited  presenta¬ 
tion  at  the  10th  International  Colloquium  on  Dynamics  of  Explosions  and  Reactive 
Systems,  Berkeley,  CA.  August,  1985. 

Droplet  Flows  with  Surface  Tension,  D.  Fyfe,  E.  Oran  and  M.  Fritts,  the  American 
Physical  Society,  Meeting  of  the  Division  of  Fluid  Dynamics,  November,  1985. 


Droplet  Deformation  and  Breakup  in  Shear  Flows,  D.E.  Fyfe,  E.S.  Oran  and  M.J.  Fritts, 
SIAM  1986  National  Meeting,  July,  1986. 


Numerical  Simulation  of  Reactive  Flows,  invited  presentation  at  NASA-Lewis  Research 
Laboratory,  Cleveland,  OH,  1986. 

Combustion  Computation,  Wright  Aeronautical  Laboratories,  invited  presentation  at 
Wright- Patterson  AFB,  OH,  1986. 

Numerical  Simulation  of  Droplets  in  Shear  Flows,  D.E.  Fyfe,  and  E.S.  Oran,  American 
Physical  Society  Meeting  of  the  Division  of  Fluid  Dynamics,  November,  1986. 

Numerical  Simulation  of  Droplet  Oscillations,  Breakup,  and  Distortion,  D.  E.  Fyfe, 
E.  S.  Oran  and  M.  J.  Fritts,  AlA.A  25th  Aerospace  Meeting,  January,  1987.  A1.\A 
Paper  No.  87-0539. 

Droplet  Deformation  and  Breakup  in  Shear  Flows,  D.E.  Fyfe,  E.S.  Oran  and  M.J.  Fritts, 
SIAM  Meeting  on  Numerical  Combustion,  March,  1987. 

Numerical  Simulation  of  Reactive  Flows,  invited  presentation  the  25th  AIA  A  Aerospace 
Sciences  Meeting,  Reno,  1987. 

Some  Challenges  in  the  Numerical  Simulation  of  Laminar  and  Turbulent  Flames,  in¬ 
vited  presentation  at  Princeton  University,  1987. 


The  book  Numerical  Simulation  of  Reactive  Flow  by  E.S.  Oran  and  J.P.  Boris,  Else¬ 
vier,  1987.  Large  portions  of  this  were  was  based  on  work  done  in  and  material 
developed  during  the  course  of  this  project. 

Surface  Tension  and  Viscosity  with  Lagrangian  Hydrodynamics  on  a  Triangular  Mesh, 
D.E.  Fyfe,  M.J.  Fritts  and  E.S.  Oran,  to  appear  in  J.  Comp.  Phys.,  1988. 

SCIENTIFIC  INTERACTIONS 

During  the  research  period,  we  worked  with  Dr.  Ann  Karagozian  from  University  of 
California,  Los  Angeles,  to  modify  her  theoretical  methods  to  describe  the  break  up 
of  droplets  due  to  shear  flows  around  the  droplet.  The  theoretical  model  includes 
surface  tension  and  viscosity  and  comes  up  with  a  droplet  break  up  criterion. 

We  have  been  talking  to  Dr.  Jack  Hansman  from  MIT,  who  has  visualizations  of  the 
breakup  of  water  droplets  in  a  shear  flow.  These  show  that  the  droplet  clearly  loses 
its  spherical  shape  before  it  goes  unstable  and  breaks  up.  His  results  are  excellent 
for  comparing  with  our  calculations. 

We  have  begun  to  interact  with  Dr.  Nassar  Ashgriz  from  the  State  L^niversity  of  New 
York  at  Buffalo.  He  has  excellent  experiments  which  we  can  model. 

We  have  been  talking  extensively  to  Professor  C.K.  Law.  from  the  University  of  Califor¬ 
nia,  Davis,  about  simulating  some  of  the  conditions  of  his  droplet-droplet  collision 
experiments. 

We  have  been  consulting  with  Dr.  Josette  Belan  from  the  Jet  Propulsion  Laboratory 
about  a  range  of  parameters  over  which  it  would  be  of  use  to  do  a  parametric 
study  of  droplet  breakup. 

During  fiscal  year  1985,  we  started  a  dialogue  with  Dr.  Murial  Ishakawa,  who  is  a 
fluid  theorist  for  one  of  the  N.4SA  shuttle  flights.  The  experiments  carried  out 
were  by  Dr.  Taylor  Wang,  of  JPL.  who  wcis  interested  in  droplet  oscillations  and 
oscillations  of  droplets  within  droplets.  Our  model  was  ideal  for  doing  calculations 
of  his  experiments,  but  due  to  lack  of  personnel  and  funding,  no  calculations  of 


their  experiments  were  ever  completed.  The  net  results  was  that  we  now  have 
some  algorithms  for  droplet  rotation  and  droplets  within  droplets,  but  they  have 
not  been  applied  to  any  problem. 

During  the  research  effort,  Gopal  Patnaik  star-ed  working  with  our  group  at  NRL.  He 
did  his  dissertation  with  Dr.  William  Sirignano,  now  at  University  of  California. 
Irvine.  Dr.  Patnaik  developed  and  did  the  calculations  for  the  droplet  model  and 
has  consulted  on  our  droplet  program. 
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Abstract 

Numerictd  algorithms  for  surface  tension  tmd  viscos¬ 
ity  are  presented  in  the  context  of  a  Lagrangian  treatment 
of  incompressible  hydrodynamics  with  a  dynamically  re¬ 
structuring  grid.  A  test  problem  consisting  of  an  oscillat¬ 
ing  droplet  is  described.  Calculations  of  droplet  distortion 
and  breakup  in  an  initially  constant  flow  and  in  an  initial 
external  shear  flow  are  presented. 

Introduction 

In  principle,  a  Lagrangian  formulation  of  the  hydro¬ 
dynamics  equations  is  particularly  attractive  for  numerical 
csdculations.  Each  discretized  fluid  element  is  tracked  as 
it  evolves  through  the  interaction  with  its  changing  envi¬ 
ronment  and  with  external  forces.  The  local  interactions 
can  be  represented  without  nonphysical  numerical  diffu¬ 
sion.  Conservation  laws  are  simple  to  express  since  there 
are  no  fluxes  out  of  the  fluid  element  boundaries.  The 
paths  of  the  fluid  elements  are  themselves  a  flow  visual¬ 
ization.  It  thus  appears  to  be  the  natural  approach  to 
transient  hydrodynamics  with  free  surfaces,  interfaces,  or 
sharp  boundaries. 

In  practice,  the  use  of  Lagrangian  methods  in  nu¬ 
merical  simulations  has  generally  been  restricted  to  “well- 
behaved”  flows.  Shear,  fluid  separation,  or  even  large  am¬ 
plitude  motion  produce  severe  grid  distortion.  These  dis¬ 
tortions  arise  because  grid  points  can  move  far  enough  that 
their  nettr- neighbors  change  in  the  course  of  a  calculation. 
When  differentisd  operators  are  approximated  over  a  mesh 
which  is  distorting,  the  approximations  may  become  inac¬ 
curate.  Attempting  to  regtun  accuracy  through  regridding 
and  interpolating  physical  quantities  onto  the  new  grid  in¬ 
troduces  numerictd  diffusion  into  the  ctilculation. 

We  first  describe  the  numerical  technique  for  La¬ 
grangian  calculations  using  a  restructuring  triangular  mesh 
[1]  for  incompressible,  two-dimensional  Cartesian  Flows. 
The  major  advance  of  this  approach  is  that  the  grid  auto¬ 
matically  tuiapts  and  refines  itself  to  maintain  accuracy  for 
discretized  operators  in  a  manner  that  is  nondiffusive.  The 
algorithms  have  been  implemented  in  the  code  SPLISH, 
which  has  been  applied  to  calculations  of  nonlinear  waves 
[2.  3],  flows  over  obstacles  |4j,  Kelvin- Helmholtz  instabili¬ 
ties  [5],  Rayleigh- Taylor  instabilities  [6j,  Couette  flows  and 
Taylor  vortex  flows  [7].  We  then  describe  the  new  surface 
tension  algorithms,  and  describe  applications  of  the  code 
to  calculations  of  droplet  oscillations  and  distortions  due 
to  the  presence  of  background  flows. 


Basic  Elements  of  Lagrangian  Trianmiiay  r.rj^c 

Consider  a  two-dimensional  space  which  is  divided 
into  triangular  cells.  A  section  of  this  mesh  shown  in 
Figure  1,  which  shows  an  interface  between  fluid  type  I 
and  fluid  type  II.  In  Figure  la.  a  particular  triangle  j  is 
highlighted  by  heavy  lines  and  the  various  components  of 
the  triangle  are  labeled.  Three  vertices,  V'l,  V2,  and  V',, 
are  connected  consecutively  by  sides  5j,  S2,  and  S3.  The 
direction  of  labeling  around  each  triangle  is  counterclock¬ 
wise  and  the  z  axis  is  directed  out  of  the  page.  Since  the 
mesh  can  be  irregularly  connected,  an  arbitrary  number 
of  triangles  can  meet  at  each  vertex.  We  can  define  a 
cell  surrounding  a  vertex,  as  shown  in  Figure  lb,  by  the 
shaded  region  surrounding  V3 .  The  borders  of  such  vertex- 
centered  cells  are  determined  by  constructing  line  segments 
joining  the  centroid  of  each  triangle  with  the  midpoints  of 
the  two  triangle  sides  connected  to  the  vertex,  for  all  tri¬ 
angles  surrounding  that  vertex.  This  definition  of  a  vertex 
cell  equally  apportions  the  area  of  a  triangle  to  each  of  its 
three  vertices  and  provides  a  simple,  efficient  way  to  evalue 
the  finite  difference  operators. 
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Figure  1.  A  section  of  a  triangular  grid  show  ing  a) 
a  material  interface,  b)  a  vertex-cell. 
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The  approach  we  describe  is  a  conservative  integral, 
control  volume  approach  on  a  triangultir  grid  that  uses  an 
integral  formulation  to  derive  the  difference  algorithms. 
We  use  the  expressions  for  the  integral  of  the  gradient  of  a 
scalar  funtion,  /,  and  the  divergence  and  curl  of  a  vector 
field,  V.  in  two  Cttrtesion  dimensions, 


/  V/  d.4  = 

(f  f  d\x  z 

(1) 

J  A 

Jc 

f  V  \dA  = 

®  V  -  ( tfl  X  i ) 

(2) 

J  A 

Jc 

/  V  X  vd.4  = 

y  ■  d\  i  , 

(3) 

'a 

Jc 

where  .4  is  the  region  enclosed  by  the  curve  C  and  d\  is 
the  vector  tire  length  around  C  in  the  counterclockwise 
direction. 

A  triangle-centered  quantity  is  assumed  to  be  piece- 
wise  constant  over  the  triangles  with  discontinuities  occur¬ 
ring  at  the  triangle  sides,  and  a  vertex-centered  quantity 
is  assumed  to  be  piecewise  linear  over  the  triangles.  If  we 
want  to  form  a  triangle-centered  derivative,  we  use  the  tri¬ 
angle  as  the  area  A  and  the  sides  of  the  triangle  for  the 
curve  C  in  Eqs.  (1)  -  (3).  If  we  want  to  form  a  vertex- 
centered  derivative,  we  use  the  vertex-centered  cell  as  the 
area  .4.  We  approximate  the  area  integral  on  the  left  side 
of  Eq.  ( 1 )  -  (3)  by  the  area  of  the  vertex-centered  cell  times 
the  value  of  the  derivative  at  the  vertex.  We  approximate 
the  line  integral  using  the  value  on  each  triangle  and  the 
appropriate  vector  length  through  the  triangle.  This  ap¬ 
proach  is  described  in  more  detail  in  Fyfe  et  aj.  [8). 

The  basic  equations  for  inviscid  incompressible  hydro¬ 
dynamics  axe 
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The  vertex  velocity  v,  in  Eq.  (11)  is  obtained  from  the 
area-weighted  Vj  determined  in  the  previous  iteration. 

The  transformation  R  results  from  the  requirement  of 
conservation  of  circulation.  Equation  ( 13),  which  produces 
conservation  of  circulation  over  vertex  cell  volumes,  is  a 
consequence  of  this  approach.  It  reflects  numerically  the 
fact  that  the  triangle  velocities  must  be  altered  as  the  grid 
rotates  and  stretches.  The  transformation  R  is  derived 
by  considering  the  circulation  about  each  vertex.  Con¬ 
servation  of  vorticity  then  takes  the  form  of  the  operator 
R  which  preserves  the  value  of  the  circulation  about  each 
vertex  as  the  grid  changes. 

The  pressures  {p"}  in  Eq.  (14)  are  derived  from 
the  condition  that  the  new  velocities  {v"}  should  be 
divergence-free  at  the  new  timestep,  satisfying  Eq.  (5). 
The  pressure  Poisson  equation  is  derived  from  Eq.  (14) 
by  setting  (V  ■  v")  =  0  to  obtain  a  pressure  p",  such  that 

(V  ■  =  (V  •  v/^),  +  (V  .  .  15) 

Both  terms  in  Eq.  (15)  are  straightforward  to  evaluate, 
since  the  divergence  is  taken  over  triangle-centered  quan¬ 
tities.  Two  features  of  the  Poisson  equation,  Eq.  (15).  are 
noteworthy.  First,  it  is  derived  from  =  V  ■  Vo,  as 
in  the  continuum  c«e.  Second,  the  left-hand-side  results 
in  the  more  familiar  second-order  accurate  templates  for 
the  Laplacitms  (such  as  the  five- point  formula)  derived  for 
homogeneous  fluids  and  regular  mesh  geometries. 


In  two  dimensions  the  fluid  density  p,  pressure  p,  and  ve¬ 
locity  V  are  assumed  to  vary  with  z,  y,  and  (.  The  term  f, 
represents  external  forces  applied  to  the  fluid,  for  extunple, 
forces  due  to  gravity. 

In  this  formulation,  it  is  important  to  consider  which 
of  the  physical  variables,  p,  v,  and  p,  should  be  defined 
as  vertex-centered  quantities  and  which  should  be  de¬ 
fined  as  triangle-centered  quantities.  Choosing  these  cor¬ 
rectly  ensure  the  correct  conservation  properties.  We  have 
found  that  prescribing  velocities  as  triangle-centered  quan¬ 
tities  makes  the  formulation  of  conservation  of  circulation 
straightforward.  Prescribing  the  densities  on  triangles  and 
pressures  at  vertices  allows  conservation  of  vertex  cell  ar¬ 
eas. 

The  numerical  integration  procedure  for  velocities 
uses  a  split-step  algorithm.  The  velocities  are  advanced 
a  half  timestep,  the  grid  is  advanced  a  full  timestep,  and 
then  the  velocities  are  advanced  forward  the  other  half 
timestep: 


Viscosity  modifies  Eq.  (6),  so  that  now 


dv 


p^  -t-  Vp  =  fe  -(-  pV*v. 


(16) 


Discretization  of  the  additional  term  in  the  momentum 
equation  follows  the  seime  approach  as  the  discretization  of 
the  other  terms.  Since  the  velocity  is  a  triangle-centered 
quantity,  we  need  a  discrete  vertex-centered  gradient  oper¬ 
ator  and  a  discrete  triangle-centered  divergence  operator. 
The  finite  difference  equations.  Eqs.  (10)  and  (14),  can  be 
modified  to  account  for  the  additional  term  in  the  momen¬ 
tum  equation  by 
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These  equations  are  implicit  in  the  velocities,  just  as  the 
original  Eqs.  (10)  -  (14)  are.  .As  in  the  inviscid  case,  we 
solve  by  iteration. 


In  Lagretngian  calculations,  the  grid  may  distort  to  the 
point  where  grid  restructuring  is  necessary.  The  deriva¬ 
tions  of  the  reconnection  and  vertex  addition  tmd  deletion 
tdgorithms  are  done  through  the  control  volume  approach 
and  the  use  of  triangle  velocities.  For  all  the  algorithms 
used,  the  divergence  and  curl  tedcen  about  each  vertex  are 
both  identically  conserved  for  grid  reconnections  and  ver¬ 
tex  addition. 

The  accuracy  of  a  general  triangular  mesh  is  dimin¬ 
ished  by  large  obtuse  angles  within  triangles.  With  recon¬ 
nections,  accuracy  can  be  recovered  by  ensuring  that  large 
obtuse  angles  are  preferentially  eliminated.  Of  the  many 
ways  of  formulating  a  reconnection  algorithm,  we  have  cho¬ 
sen  one  based  on  requirements  for  solving  the  pressure 
Poisson  equation.  Since  the  equation  is  solved  by  itera¬ 
tion,  we  want  the  iteration  to  converge  as  rapidly  as  pos¬ 
sible.  Mathematictdly,  convergence  is  assured  if  the  finite 
difference  equation  is  diagonally  domineuit.  This  require¬ 
ment  translates  to  a  relation  between  two  of  the  angles  of 
each  triangle.  The  reconnection  algorithm  preferentially 
eliminates  large  angles  in  triangles,  since  the  diagontd  is 
chosen  which  divides  the  largest  opposing  angles.  Inter¬ 
face  sides  are  never  allowed  to  reconnect.  In  such  cases 
vertex  addition  algorithms  are  needed. 

Vertex  addition  algorithms  axe  tdso  needed  where  the 
flow  naturally  depletes  vertices.  For  vertex  addition,  satis¬ 
faction  of  conservation  integrals  is  particularly  simple.  The 
vertex  added  at  the  centroid  of  a  triangle  subdivides  that 
triangle  into  three  smaller  triangles.  A  vertex  added  to  the 
midpoint  of  a  side  subdivides  the  two  adjacent  triangles 
into  four  smaller  triangles.  If  the  new  triangle  velocities 
are  all  the  same  as  the  velocity  of  the  subdivided  triangles, 
all  conservation  laws  are  satisfied.  Since  the  reconnection 
algorithm  is  also  conservative,  subsequent  reconnections  to 
other  vertices  ensure  that  the  only  effect  of  the  addition  is 
an  increase  in  resolution. 

The  case  is  not  as  obvious  for  vertex  deletion.  Recon¬ 
nections  can  be  used  to  surround  any  interior  vertex  within 
a  tritmgle.  The  vertex  is  then  removed  and  the  new  lEu-ger 
triangle  given  a  velocity  which  is  the  area-weighted  sum  of 
the  old  velocities,  redistributes  circulation  in  accordance 
with  area  coordinates. 


The  surface  tension  at  an  interface  between  two  ma¬ 
terials  depends  on  the  curvature  of  the  interface.  In  the 
conventional  numerical  representation  of  surface  tension, 
it  is  cast  into  a  finite-difference  form  by  fitting  vertices  on 
the  material  interface  to  some  parametric  function.  This 
function  is  then  used  to  find  an  estimate  of  local  curva¬ 
ture.  Once  the  curvature  is  known,  a  surface  tension  force 
is  evaluated  and  used  to  accelerate  interface  vertices. 

This  scheme  fails  in  SPLISH  for  two  reasons.  First, 
the  interface  vertices  are  accelerated  directly  by  surface 
tension  forces  evaluated  on  the  vertices.  Since  velocities 
au’e  centered  on  triangles  in  SPLISH,  the  velocity  field  sees 
the  effect  of  the  acceleration  a  half-timestep  later,  unless 
a  secondary  calculation  is  made.  As  a  result,  the  pres¬ 
sure  calculated  within  the  droplet  is  inconsistent  with  that 


found  from  the  surface  tension  formula.  Second,  since 
the  pressure  gradient  forces  and  surface  tension  forces  are 
not  calculated  in  the  same  manner,  numerical  errors  result 
which  grow  with  each  timestep. 

Both  of  these  problems  are  eliminated  by  a  different 
formulation  of  surface  tension,  in  which  a  surface  tension 
potential  is  used  to  generate  the  forces.  The  surface  tension 
force  is  formulated  as  a  gradient  of  a  potential  present  only 
at  the  surfaces.  With  this  method,  the  pressure  gradient 
forces  are  calculated  in  the  same  manner  tmd  on  the  same 
grid  ais  the  forces  derived  from  the  surface  tension  poten¬ 
tial.  Therefore  both  the  surface  tension  potential  and  the 
pressure  are  dynamically  similar,  and  the  physical  pressure 
drop  across  the  interface  must  exactly  cancel  the  surface 
tension  forces.  Preliminary  aspects  of  this  work  were  de¬ 
scribed  by  Fritts  et  al.  [8.9], 

The  finite-difference  algorithms  for  surface  tension  are 
straightforward.  The  surface  tension  forces  are  included 
through  Laplace’s  formula  for  the  pressure  jump  across  an 
interface  (llj, 

p,-p„=a/R.  (16) 

where  p,  is  the  pressure  just  inside  the  droplet  at  the  inter¬ 
face,  Pa  is  the  pressure  just  outside  the  droplet  at  the  inter¬ 
face,  a  is  the  surface  tension  coefficient  sissociated  with  the 
two  media  which  define  the  interface,  and  R  is  the  radius 
of  curvature  in  the  two-dimensional  plane.  The  radius  of 
curvature  is  positive  at  points  on  the  interface  where  the 
droplet  surface  is  convex  (a  circle  is  convex  everywhere) 
and  negative  when  the  surface  is  concave.  These  pressure 
jumps  are  included  in  the  Poisson  equation  for  the  pres¬ 
sure.  The  average  pressure.  (Pi  -(-Po)/2.  is  computed  at  the 
interface  vertices.  From  the  average  pressure  and  the  pres¬ 
sure  jump,  we  can  compute  a  pressure  gradient  centered 
on  triangles,  both  inside  and  outside  the  surface.  This 
pressure  gradient  is  used  in  the  momentum  equation. 

The  radius  of  curvature  is  computed  from  a  paramet¬ 
ric  cubic  spline  interpolant  to  the  interface  vertices.  Past 
calculations  of  droplets  oscillating  due  to  surface  tension 
forces  (12,13)  also  use  cubic  spline  interpolation.  However, 
they  divided  the  surface  into  at  least  four  segments  (the 
top,  bottom,  right  and  left  sides  of  the  droplet)  to  pro¬ 
duce  an  interpolant  on  each  segment.  Each  interpolant 
was  matched  at  the  joints  to  produce  an  overall  curve. 
The  parametric  interpolant  used  here  does  not  require  this 
special  matching.  We  generate  the  twice  differentiable  pe¬ 
riodic  spline  interpolants.  r(s)  =  (r(s),  j/ls  l)  as  prescribed 
by  deBoor  (14). 

The  spline  fit  is  also  used  for  regridding.  When  the 
regridding  algorithm  calls  for  the  bisection  of  a  triangle 
side  which  borders  the  two  media,  a  new  vertex  is  added 
on  the  spline  interpolant  between  the  vertices.  This  is 
done  rather  than  bisecting  the  straight-line  segment,  since 
a  straight-line  bisection  introduces  spurious  interface  os¬ 
cillations.  Bisecting  the  spline  meiinttuns  a  better  overall 
shape  for  the  interface. 
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In  order  to  test  the  2dgorithm  for  surface  tension  in 
SPLISH.  we  calculated  the  oscillation  of  a  droplet  due  to 
surface  tension.  We  have  extended  Rayleigh's  linear  theory 
[15]  of  small  amplitude  droplet  oscillations  to  include  the 
presence  of  an  external  fluid. 


=  (n’  —  n) 


iPd  +  P«)a^  ’ 


where-.'  is  the  frequency,  is  the  droplet  density,  p,  is  the 
density  of  the  external  fluid,  and  the  surface  of  the  droplet 
is  given  in  polar  coordinates  by 

r  =  a  +  ecos(nS),  (18) 

where  a  is  the  unperturbed  radius  of  the  jet,  and  n  pre¬ 
scribes  the  mode  of  oscillation  in  the  plane  with  amplitude 


Figure  2  is  a  composite  of  frames  from  a  calculation  in 
which  e  =  0.2a  =  0.0025.  In  this  calculation  there  are  17 
vertices  in  each  direction  along  the  exterior  boundaries.  12 
vertices  on  the  droplet  interface  and  a  total  of  313  vertices 
initirdly  in  the  calculation.  The  computational  domain  is 
0.1  cm  on  a  side.  The  left  cind  right  boundaries  are  peri¬ 
odic  while  the  top  and  bottom  boundaries  tire  solid  walls. 
The  timestep  is  dt  =  2.5  x  10“^  s.  The  figures  shows  four 
and  a  half  oscillations  of  the  droplet.  We  can  see  that  as 
the  calculation  proceeds,  no  new  vertices  have  been  added, 
but  in  fact  some  have  been  subtracted.  This  was  the  case 
because  the  initial  gridding  was  adequate  to  represent  the 
droplet  shape.  From  these  ctilculations,  the  period  of  os¬ 
cillation  is  r,2  =  1.35  X  10“'’  s,  compared  to  the  theoretical 
value  of  1.13  x  10”’  s.  Figures  3  shows  the  initial  oscilla¬ 
tion  for  a  more  resolved  case  in  which  there  tire  28  vertices 
surrounding  the  droplet. 


Incompressible  Flow  about  a  Dropipt 

In  this  section  we  present  some  calculations  of  forced, 
asymmetric  drop  oscillations  induced  by  flow  around  a 
droplet.  These  calculations  include  both  the  effects  of  vis¬ 
cosity  and  surface  tension.  The  capability  of  studying  such 
flows  for  viscous  droplets  in  shear  flows  is  the  motivation 
for  developing  the  viscosity  and  surface  tension  algorithms. 

The  initial  conditions  we  used  specified  an  initially 
steady-state  potential  flow  about  a  periodic  series  of  cylin¬ 
ders.  .\gain,  the  boundary  conditions  on  the  left  and  right 
sides  are  periodic,  and  the  upper  and  lower  boundary  con¬ 
ditions  are  reflecting  walls.  Initially,  a  perfectly  circular 
droplet  is  at  rest  in  a  background  flow.  A  physical  situa¬ 
tion  modelled  by  such  an  initialization  might  occur  if  the 
flow  velocity  were  ramped  up  to  its  final  value  before  anv 
significant  structure  could  develop  in  the  flow,  and  before 
the  droplet  could  pick  up  any  substantial  velocity.  Basi¬ 
cally,  it  is  a  smooth  stMt  for  the  calculation.  Previouslv  we 
had  performed  ctJculations  which  began  with  an  impulsive 
start,  but  found  that  as  a  result  there  was  a  large  amount 
of  momentum  transferred  across  the  droplet  interface  earlv 
in  the  calculation. 


Table  1.  Conditions  for  Flow  around  Droplet  Calculation 


density  of  kerosene 
density  of  air 
surface  tension  (STPl 
viscosity  of  kerosene 
viscosity  of  air 
air  velocity 

initial  droplet  velocity- 
droplet  radius 


0.S2  g/cc 
0.0013  g/cc 
30  dynes /cm 
1.8  centipoise 
0.018  centipoise 
120  m/s 
0.0  m/s 
125  microns 


Figure  4.  Pathlines  from  a  calculation  of  a  flow  around  a  kerosene  droplet  at  flow  velocity 
of  120  m/s  and  Re  a:  2000. 


k-»  ii->  ■/>  i.- .  J  wjrt^  »>'  »Jf  J  ’>  'y  >.’■"  J  '■  Ji  ■  JJ  V  V  *j 


The  cedculations  presented  here  model  the  forced  fluid 
flow  due  to  a  fast  air  stream  about  an  initially  stationary 
kerosene  droplet.  The  physical  parameters,  given  in  Ta¬ 
ble  1,  are  appropriate  for  a  combustor  environment.  A 
total  of  309  vertices  were  used  to  initialize  the  problem, 
w*  with  12  vertices  at  the  droplet  interface.  Figure  4  follows 
the  evolution  of  pathlines  in  the  internal  and  external  flow 
fields  through  a  series  of  timesteps.  For  an  air  velocity  of 
120  m/s  and  a  droplet  radius  of  125  microns,  the  corre¬ 
sponding  Reynolds  number  is  roughly  2000.  The  pathlines 
are  defined  by  the  paths  of  vertices  over  five  timesteps. 
By  the  last  frame  of  Figure  4.  the  fluid  originally  to  the 
left  of  the  droplet  htis  progressed  through  the  mesh  and 
interacted  with  the  face  of  the  (next)  droplet. 

The  first  clear  indication  of  the  development  of  the 
recirculation  region  is  seen  in  the  fourth  frame  of  Figure  4. 
which  shows  a  pair  of  counter-rotating  vortices.  The  recir¬ 
culation  zone  continues  to  develop  throughout  the  calcu¬ 
lation,  although  at  times  the  vortex  pair  is  not  eis  evident 
due  to  the  deletion  and  addition  of  vertices,  which  inter¬ 
rupts  the  continuity  of  the  pathlines.  By  the  last  frame, 
another  pair  of  vortices  is  forming  near  the  droplet,  and 
the  original  ptdr  has  been  shed. 


Distortions  in  the  face  of  the  droplet  are  evident  in  at 
least  the  seventh  frame.  These  distortions  occur  because 
the  curvature  has  increased  and  the  streamlines  m  the  ex¬ 
ternal  flow  aie  condensed  by  the  approaching  wake.  The 
internal  velocities  are  smtdl  compared  to  the  externtil  flow 
rates  and  therefore  cannot  be  distinguished  as  pathlines 
However,  indication  of  the  (small)  internal  reciriulation 
may  be  obtaiined  by  comparing  internal  vertex  position.^ 
at  various  timeteps. 

Figure  5  shows  the  grid  at  times  in  the  calculation 
corresponding  to  those  in  Figure  4.  During  the  course  of 
the  calculation,  a  great  deed  of  vertex  addition  and  deletion 
has  occured.  Vertex  addition,  however,  is  most  noticeable 
in  the  wake  of  the  droplet  and  around  the  droplet  interface. 
Whereas  there  were  300  vertices  at  the  beginning  of  the 
calculation,  there  are  450  at  the  end. 

As  seen  in  Figure  5.  the  computational  grid  needs  fur¬ 
ther  refinement  at  this  time  because  the  perturbations  can¬ 
not  be  resolved  by  the  limits  set  on  minimum  triangle  size 
originzdly  chosen  for  the  ceJculation.  A  sign  that  the  cal¬ 
culation  IS  under- resolved  is  that  one  of  the  crests  of  the 
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Figure  5.  Frames  showing  the  triangular  grid  at  the  same  times  as  shown  for  the  pathlines 
in  Figure  4. 
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surface  wave  is  spanned  by  a  single  triangle,  a  situation 
which  allows  no  communication  of  that  surface  fluid  with 
the  interior  of  t he  droplet .  In  order  to  continue  the  simula- 
tion,  better  resolution  must  be  obtained  about  the  droplet 
surface.  Another  algorithm  is  currently  being  included  to 
allow  higher  resolution  near  points  of  large  curvature  at 
^  material  interfaces. 


■A  Droplet  in  a  Shear  Flow 


An  important  problem  in  atomization  is  how  a  droplet 
breaks  up  due  to  shear  forces.  To  investigate  this  compu¬ 
tationally,  we  have  simulated  droplets  in  a  shear  flow  cen¬ 
tered  around  the  droplet.  The  initial  flow  is  prescibed  bv 
t’l  =  G(y  —  yd)  for  points  outside  the  droplet  and  c,  =  0 
for  points  inside  the  droplet.  The  parameter  yd  is  the  y- 
coordinate  of  the  center  of  the  drop  and  G  gives  the  mag¬ 
nitude  of  the  shear.  The  y-component  of  the  velocity  is 
initially  zero  everywhere. 


Results  of  a  such  a  shear  on  a  kerosene  droplet  in  hot 
air  axe  shown  in  the  three  panels  in  Figures  6.  for  a  case  in 
which  G  =  3  X  10^/sec.  We  used  0.013y/cc  as  the  density 
for  hot  air.  The  remaining  physical  parameters  are  the 
same  as  those  in  Table  1.  Initially  the  droplet  was  round, 
but  in  Figure  6a  it  has  already  become  elongated  in  the 
directions  of  the  shear.  .At  a  later  time,  in  Figure  6b,  it 
has  become  even  more  elongated.  Times  intermediate  to 
these  two  figures  show  that  some  very  small  droplets  have 
already  been  pulled  off  of  the  large  drop,  but  their  size  wais 
so  small  that  they  were  deleted  from  the  calculations.  Fig¬ 
ure  6c  shows  a  still  later  time,  when  smtill  droplets  have 
been  are  moving  off  of  both  sides.  The  small  droplets  some¬ 
times  seem  to  move  counter  to  the  flow  of  the  main  shear 
layer.  This  is  because  a  recirculation  zone  forms  on  the 
upper  left  and  the  lower  right  of  the  large  droplet. 
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Summary 

This  paper  presented  the  current  algorithm.s  included 
m  the  code  SPLISH.  a  two-dimensional  Cartesian  La- 
grangian  treatment  of  incompressible  flows  with  a  dynam¬ 
ically  restructuring  grid.  .Algorithms  for  modelling  viscos¬ 
ity  and  surface  tension  have  been  tested  on  a  number  of 
problems. 

Detailed  benchmarks  of  the  surface  tension  algorithm 
were  presented  using  a  Rayleigh  oscillating  droplet  for  a 
test  problem.  This  algorithm,  based  on  spline  fits  to  de¬ 
termine  curvature,  was  good  enough  to  allow  the  droplet 
to  oscillate  mtmy  times  and  still  maintain  a  constant  pe¬ 
riod.  However,  the  amplitude  calculated  for  the  original 
excited  mode  decayed  into  higher  modes.  We  believe  this 
is  due  to  a  resonance  coupling  between  the  modes  n  =  2 
and  n  =  3.  .A  similar  resonance  coupling  exists  in  three 
dimensions[16j. 

We  au-e  currently  calculating  droplet  distortion  and 
breakup  due  to  differential  externad  flows,  shear  flows,  and 
droplet-droplet  collisions.  Some  of  the  results  of  these  cal¬ 
culations  were  shown  m  Figures  4.  5.  and  6.  and  will  be 
discussed  more  thoroughly  in  the  presentation. 
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Abstract 


Numerical  algorithms  for  surface  tension  and  viscosity  are  presented  in  the  con¬ 
text  of  a  Lagrangian  treatment  of  incompressible  hydrodynamics  with  a  dynamically 
restructuring  grid.  New  algorithms  are  given  which  update  previous  Lagrangian  ap¬ 
proaches  in  the  code  .SPLISH.  Test  problems  involving  internal  gravity  and  capillary 
waves,  an  oscillating  droplet  and  a  viscous  shear  layer  are  described,  .^n  example  is 
given  of  a  flow  calculated  in  and  around  a  viscous  droplet  with  surface  tension  in  a 
shear  flow. 


>v 


I.  Introduction 


0-5C.' 


A 


In  principle,  a  Lagrangian  formulation  of  the  hydrodynamics  equations  is  partic¬ 
ularly  attractive  for  numerical  calculations.  Each  discretized  fluid  element  is  tracked 
as  it  evolves  through  the  interaction  with  its  changing  environment  and  with  e.xternal 
forces.  The  local  interactions  can  be  represented  without  nonphysical  numerical  diffu¬ 
sion.  Conservation  laws  are  simple  to  e.xpress  since  there  are  no  flu.xes  out  of  the  fluid 
element  boundaries.  The  paths  of  the  fluid  elements  are  themselves  a  flow  visualiza¬ 
tion.  It  thus  appears  to  be  the  natural  approach  to  transient  hydrodynamics  with  free 
surfaces,  interfaces,  or  sharp  boundaries. 

In  practice,  the  use  of  Lagrangian  methods  in  numerical  simulations  has  generally 
been  restricted  to  "well-behaved”  flows.  Shear,  fluid  separation,  or  even  large  amplitude 
motion  produce  severe  grid  distortion.  These  distortions  arise  because  grid  points  can 
move  far  enough  that  their  near-neighbors  change  in  the  course  of  a  calculation.  When 
differential  operators  are  approximated  over  a  mesh  which  is  distorting,  the  approx¬ 
imations  may  become  inaccurate.  .Attempting  to  regain  accuracy  through  regridding 
and  interpolating  physical  quantities  onto  the  new  grid  introduces  numerical  diffusion 
into  the  calculation. 

This  paper  is  a  summary  and  update  of  the  latest  additions  and  modifications 
to  a  numerical  technique  for  indefinitely  extending  Lagrangian  calculations  by  using 
a  restructuring  triangular  mesh,  first  introduced  by  Fritts  and  Boris  [1].  The  major 
advance  of  this  approach  is  that  the  grid  automatically  adapts  and  refines  itself  to  main¬ 
tain  accuracy  for  discretized  operators  in  a  manner  that  is  nondiffusive.  The  algorithms 
have  been  implemented  in  the  code  SPLISII.  which  has  l)een  applied  to  ralcnlal  ions 
of  nonlinear  waves  [2.  3],  flows  over  obstacles  [-1],  Kelvin-llehnholt z  inst abilit  it's  i")|. 
Rayleigh-Taylor  instabilities  [6].  Couette  flows  and  Taylor  vortex  flows  [7]. 

Work  on  Lagrangian  techniques  for  grids  which  do  not  have  fixed  connectivii \- 
has  recently  had  a  renaissance.  Early  attempts  includetl  the  P.\X.\CE.\  code  [Si  and 
the  P.AF  ( Particle-.And-Force)  algorithm  [9.  10].  In  the  I'JTO's.  these  concepts  wer<' 
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improved  and  extended  for  triangular  grids:  triangle  reconnection  Ijv  Crowley  [11]: 
MHD  algorithms  over  a  triangular  mesh  [12];  and  adaptive  triangular  meshes  in  the 
work  mentioned  in  the  previous  paragraph  on  SPLISH,  During  the  same  period  work 
began  which  used  Voronoi  meshes  for  hydrodynamics  calculations  [13|. 

Recently  this  use  of  general  connectivity  grids  has  rapidly  expanded,  as  sumnta- 
rized  in  the  First  International  Conference  on  Free-Lagrange  Methods  [14|.  Applica¬ 
tions  now  include  finite-difference  and  finite-element  calculations  of  classic  hydrody¬ 
namic  instabilities,  tokamak  modelling,  high  temperature  plasma  physics,  heat  con¬ 
duction,  wave-structure  interactions,  impact  deformations,  and  hydrodynamics  prob¬ 
lems  for  both  compressible  and  incompressible  fluids.  Free-Lagrange  methods  now  use 
quadrilateral,  triangular  and  mixed  meshes  in  two  dimensions,  tetrahedral  meshes  in 
three  dimensions,  Moronoi  meshes  in  both  two  and  three  dimensions,  and  methods 
which  are  mesh-free. 

In  this  paper  we  present  the  latest  modifications  to  SPLISH  (section  II).  These 
include  the  most  recent  version  of  the  rotation  operator,  which  conserves  circulation, 
and  the  residual  algorithm,  which  ensures  conservation  of  the  area  of  cells.  We  also 
introduce  new  algorithms  for  viscosity  and  surface  tension.  Including  viscosity  proved 
to  be  straightforward  (section  II).  However,  the  search  for  a  good  enough  algorithm 
for  surface  tension  (section  III)  was  more  challenging  and  difficult.  The  basic  problem 
is  defining  a  proper  curvature  from  a  finite  number  of  points.  Because  of  this,  the 
numerical  approximation  of  surface  tension  forces  between  two  fluids  is  conceptually 
quite  different  from  approximations  of  convection  and  viscous  forces.  The  final  formu¬ 
lation  chosen,  a  series  of  test  problems,  and  a  list  of  approaches  that  failed  are  detailed 
(section  III).  Finally,  we  combine  the  convective  transport,  surface  tension,  and  vis¬ 
cosity  algorithms  to  perform  some  preliminary  calculations  of  flows  in  and  around  a 
viscous  kerosene  droplet.  These  calculations  show  vortex  shedding  behind  the  droplet, 
distortion  of  tlie  droplet  due  to  the  shear  flow,  and  internal  droplet  flows. 


II.  Basic  Elements  of  Lagrangian  Triangular  Grids 

This  section  is  a  review  of  the  derivation  of  low  order  finite-difference  approxi¬ 
mations  to  the  equations  describing  incompressible  fluid  motion  for  general  triangular 
grids.  Some  of  the  material  was  originally  presented  l)y  Fritts  and  Boris  [1],  and  the 
interested  reader  is  referred  there  for  more  detail.  However,  new  material  brings  the 
previous  paper  up  to  date.  This  includes  the  lastest  version  of  the  rotation  operator, 
which  conserves  circulation,  the  residual  algorithm,  which  ensures  conservation  of  the 
area  of  cells,  and  the  new  algorithm  for  viscosity. 

.-4.  The  Triangular  Grid 

Consider  a  two-dimensional  space  which  is  divided  into  triangular  cells.  A  section 
of  this  mesh  shown  in  Fig.  1.  which  shows  an  interface  between  fluid  type  I  and  fluid 
type  II.  In  Fig.  la,  a  particular  triangle  j  is  highlighted  by  heavy  lines  and  the  various 
components  of  the  triangle  are  labeled.  Three  vertices.  Vj,  lA,  and  V3,  are  connected 
consecutively  by  sides  5i,  S2,  and  S3.  The  direction  of  labeling  around  each  triangle 
is  counterclockwise  and  the  ;  axis  is  directed  out  of  the  page.  Since  the  mesh  can  be 
irregularly  connected,  an  arbitrary  number  of  triangles  can  meet  at  each  vertex. 

We  can  define  a  cell  surrounding  a  vertex,  as  shown  in  Fig.  lb,  by  the  shaded 
region  surrounding  Is.  The  borders  of  such  vertex-centered  cells  are  determined  by 
constructing  line  segments  joining  the  centroid  of  each  triangle  witli  the  midpoints  of 
the  two  triangle  sides  connected  to  the  vertex,  for  all  triangles  surrounding  that  vertex. 
This  definition  of  a  vertex  cell  equally  apportions  the  area  of  a  triangle  to  each  of 
its  three  vertices  and  provides  a  simple,  efficient  way  to  evalue  the  finite  difference 
operators.  However,  the  definition  of  a  vertex  cell  is  arbitrary.  Other  definitions  could 
be  equally  well  employed,  although  they  generally  require  additional  calculations  to 
determine  cell  intersection  points.  The  integration  of  cell  quantities  may  therefore 
involve  more  arithmetic  operations  for  other  definitions. 


B.  Finite  Differences  on  a  Triangular  Grid 


Finite-difference  approximations  for  derivatives  of  functions  defined  on  the  trian¬ 
gular  grid  are  derived  from  the  expressions  for  the  integral  of  the  gradient  of  a  scalar 
funtion,  /,  and  the  divergence  and  curl  of  a  vector  field,  v,  in  two  Cartesion  dimensions. 


(2.1) 

(2.2) 

(2.3) 


In  each  of  these  expressions,  A  is  the  region  enclosed  by  the  curve  C  and  dl  is  the 
vector  arc  length  around  C  in  the  counterclockwise  direction.  The  variable  i  is  a  unit 
vector  in  the  direction  of  the  ignorable  coordinate.  By  using  these  definitions  in  a 
conservative  integral  approach,  the  definitions  for  spatial  derivatives  described  below 
can  be  naturally  extended  to  two-dimensional  axisymmetric  geometry  [7]. 

Throughout  the  following  discussion  a  triangle-centered  cjuantity  is  assumed  to  be 
piecewise  constant  over  the  triangles  with  discontinuities  occurring  at  the  triangle  sides 
and  a  vertex-centered  quantity  is  assumed  to  be  piecewise  linear  over  the  triangles.  If 
we  want  to  form  a  triangle-centered  derivative,  we  use  the  triangle  as  the  area  .1  and 
the  sides  of  the  triangle  for  the  curve  C  in  Eqs.  (2.1)-(2.3).  We  then  approximate 
the  area  integral  by  the  area  of  the  triangle  times  the  value  of  the  derivative  on  the 
triangle,  and  approximate  the  line  integral  using  the  trapezoidal  rule  on  each  side  of 
the  triangle.  For  example,  the  gradient  of  a  scalar  function  /  defined  at  the  vertices  is 
a  triangle-centered  quantity,  (V/)^.  given  by 


(r.-i-r,  +  i) 


where  r,  =  (.r,.//,)  is  a  vector  coordinate  for  vertex  i  and  .1^  is  the  area  id'  triangle  j. 
We  have  also  used  the  notation  of  Fritts  and  Boris  [1]  that 


is  interpri'ted  as  i  he 


sum  over  vertices  i  of  triangle  j.  In  the  material  presented  below,  the  index  i  designates 
vertex-centered  quantities  and  the  index  j  designates  triangle-centered  quantities. 

If  we  want  to  form  a  vertex-centered  derivative,  we  use  the  vertex-centered  cell 
as  the  area  .4.  We  approximate  the  area  integral  on  the  left  side  of  Eq.  (2.1)-(2.3) 
by  the  area  of  the  vertex-centered  cell  times  the  value  of  the  derivative  at  the  vertex. 
We  approximate  the  line  integral  using  the  value  on  each  triangle  and  the  appropriate 
vector  length  through  the  triangle.  For  example,  the  curl  of  the  vector  field  v  at  a 
vertex  c  is  approximated  by 

.4e(V  X  v)c  =  -  ^  v.+i/o  ■  (r,+,  -  r.)  f.  (2.5) 

•(C) 

where  .4c  =  j  Hj(c)  ‘‘^j  the  vertex-centered  cell  area,  ^^(c)  's  a  sum  over  the  triangles 
around  the  central  vertex  c,  ^,(c)  's  a  sum  over  the  vertices  around  vertex  c.  and 
is  the  value  of  the  vector  field  v  on  the  triangle  having  vertices  c.  i.  i  -f  1.  Similarly, 
the  divergence  of  the  vector  field  v  at  a  vertex  is  approximated  by 

.4c(V  •  v)c  =  ^  ^  (■''  +  1  “  ‘b)!  •  (2.G) 

•  (C) 


C.  The  Equations  for  [ncompressible,  Inviscid  Flow 

The  basic  equations  for  inviscid  incompressible  hydrodynamics  are 


^  =  0, 
dt 

V  •  V  =  0. 
dv 

A.^  +  v„=r.. 


(2.b) 


In  two  dimensions  the  fluid  density  />.  pressure  p.  and  velocity  v  are  assumed  to  \ai\ 
with  .r.  y.  and  t.  The  term  fc  represents  external  forces  applic'd  to  the  llnid.  tor  example. 


forces  due  to  gravity.  Equation  (2.8),  the  condition  for  incompressibility,  removes  the 
sound  waves. 


Since  we  want  our  finite  difference  approximation  to  preserve  the  conservation 
properties  for  incompressible  inviscid  fluids,  it  is  important  to  consider  which  of  the 
physical  variables,  p,  v,  and  p.  should  be  defined  as  vertex-centered  quantities  and 
which  should  be  defined  as  triangle-centered  quantities.  We  have  found  that  prescrib¬ 
ing  velocities  as  triangle-centered  quantities  makes  the  formulation  of  conservation  of 
circulation  straightforward.  Prescribing  the  densities  on  triangles  and  pressures  at 
vertices  allows  conservation  of  vertex  cell  areas. 

The  time  integration  of  velocities  uses  a  second-order  implicit  split-step  algorithm 
which  is  solved  by  iteration.  The  vertex  positions  are  advanced  using  a  second-order 
midpoint  rule.  Specifically,  the  velocities  are  advanced  a  half  timestep.  the  grid  is 
advanced  a  full  timestep.  and  then  the  velocities  are  advanced  forward  the  other  half 
timestep.  The  complete  algorithm  is  as  follows.  First  compute  the  half-timestep  trian¬ 
gle  velocities  using 

v/*  =  v;  -  — (Vp);  +  — fe.  (2.10) 

’  ^  ip,  ■'  ip, 


where  the  superscript  o  designates  the  values  at  the  old  time  step.  We  then  make  an 
initial  guess  for  the  new  triangle  velocities 


\/2.k 


n  ,k 


.1/2./L- 

V 

J 


n  ,k 


(2.11) 

=  x^  +  6tvy--\ 

(2.12) 

=  R({xr}.{x;'-^})-v/T 

(2.12) 

=  v/  -—(Vp)  +— r. 

■'  ip,  •'  ip. 

(2.M) 
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where  the  second  superscript  indicates  the  iteration  number.  The  vertex  velocity  v"'^ 
in  Eq.  (2.11)  is  obtained  from  a  weighted  average  of  the  triangle  velocities  vj'^  for 
those  triangles  having  i  as  a  vertex, 


V.  = 


:2.i5) 


Wi 


We  use  Wj  —  OjpjAj,  where  6j  is  the  angle  (in  radians)  of  triangle  j  at  vertex  i  divided 
by  TT.  The  transformation  R  in  Eq.  (2.13)  results  from  the  requirement  of  conservation 
of  circulation,  and  is  discussed  in  Section  D  below. 

The  pressures  in  Eq.  (2.14)  are  derived  from  the  condition  that  the  new 

velocities  {v"'*^}  should  be  divergence-free  at  the  new  timestep,  satisfying  Eq.  (2.8). 
The  pressure  Poisson  equation  is  derived  from  Eq.  (2.14)  by  setting  (V  •  v^)"'*'  =  0  to 


obtain  a  pressure  p"’*,  such  that 


St 


(V  .  ^(Vp)^.)"’'-  =  (V-*=  •  vj/-’*).  -h  (V 


''n.k 


St 


2p, 


2pj 
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Doth  terms  in  Eq.  (2.16)  are  straightforward  to  evaluate,  since  the  divergence  is  taken 
over  triangle-centered  quantities.  Note  also  that  the  discrete  gradient  operator  V  must 
also  carry  time  advancement  superscripts  since  it  depends  on  the  current  grid  location. 
(See  Eq.  (2.4).)  Two  features  of  the  Poisson  equation,  Eq.  (2.16).  are  noteworthy.  First, 
it  is  derived  from  V-4>  =  V  •  Vc>.  as  in  the  continuum  case.  Second,  the  left-hand-side 
results  in  the  more  familiar  second-order  accurate  templates  for  the  Laplacians  (such 
as  the  five-point  formula)  derived  for  homogeneous  fluids  atid  regular  mesh  geometries. 


D.  Conservation  of  Circulation 


The  approach  we  have  outlined  is  basically  a  control  volume  approach  which  uses 
an  integral  formulation  to  derive  the  difference  algorithms.  Equation  (2.13).  which 
produces  conser\’ation  of  circulation  over  vertex  cell  volumes,  is  a  couseqticnce  of  this 
approach.  It  reflects  numerically  the  fact  tliat  the  triangle  velocities  must  le’  alteie<l 
as  the  grid  rotates  and  stretches.  This  process  does  not  prevent  the  addition  or  loss  of 
vorticity  due  to  external  forces  or  changes  in  density  at  interfaces.  Rather  it  corrects 
any  numerical  errors  that  may  arise  because  the  grid  has  moved.  Thus  it  guarantees 
conservation  of  circulation  at  those  vertices  where  the  circulation  theorem  applies. 

The  transformation  R  is  derived  by  considering  the  circulation  about  each  vertex. 
Since  triangle  velocities  are  constant  over  the  triangle,  the  circulation  taken  aboitt 
the  boundary  of  the  vertex  cell  can  be  calculated  from  Eq.  (2.5).  The  conservation 
of  vorticity  then  takes  the  forui  of  the  operator  R  which  preserves  the  value  of  the 
circulation  about  each  vertex  as  the  grid  changes. 

Conservation  of  circulation  requires  that  at  each  timestep.  and  for  each  vertex,  c. 


E 

l(c) 


-1/2. k 

V  ' 

'^i+1/2 


-  r"’* 

T,  +  i  r,- 


)  =  E 

•■(c) 


r'/-  .  fr° 

1  +  1/2  +  l 


-rD- 


(2.r 


For  convenience  in  notation,  we  now  drop  the  superscript  1/2  for  the  velocities  and 
the  iteration  superscript  k  appearing  in  Eq.  (2.17).  Since  there  are  two  componenis 
of  velocity  on  each  triangle,  but  only  one  constraint  at  each  vertex,  the  form  of  the 
rotator  is  underdetermined.  Fritts  and  Boris  [1]  provided  the  additional  constraints  by 
making  each  term  in  ihe  circulation  integral  associated  with  a  gi\'en  triangle  a  cunseiued 
(piantity,  and  hence  the  sum  in  Eq.  (2.17)  remains  unchanged.  This  means  that  fur 
each  triangle  j. 


T.  +  i  -  r.  )  = 


T-.  +  i  -  r. 


I  =  1.2.3. 


l2.1si 


.-Mthough  this  apiuuach  conserves  circulation,  the  following  exanqde  shows  ihai  it 
much  too  restricti\’c. 
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Let  us  also  decompose  ^v,_).i/2  component,  /,  +  i/j.  luuallol  to  tlie  sidi'  op¬ 

posite  vertex  c,  and  a  component.  n,+i/2-  normal  to  the  side  opposite  vertex  e  hy 


writing 


^  ;  X  (r,+i  -  r.)  ^  ^  r,  +  ,  -  r, 

«V,  +  i/2  =  tt,+  l/2— 7-  Ti  +  ''+1/- Ui  r 

|c,+  i  —  r,|  ic,  +  i  —  r, 


(2.21) 


With  this  notation  and  using  the  equation  for  the  area.  of  triangle  i  +  1  /2. 


2-4.+ 1/2  =  J  •  [(r.+i  -  r.)  X  (r"  -  r.+i)]. 


Eq.  (2.20)  becomes 


2^.+i/2  ,  (r.+i  -  r.)  •  (r"  -  r.+i)^ 

) - r^.  +  l/2  H - i - 1 - L  +  1/2 

|r.+i-r.|  |r,+i  -  r.l 


,  -2.4._i/2  ,  (r,  -  r._])  ■  (r._]  -  r'‘), 

+  1  r'Ji-l/2  T - i  i  0-1/2 

|r. -r._i|  '  |r, -r._i|  ' 


(2.23) 


=  (v,_i/2  -  v,+i/2) -ifirc. 

Let  A’c  denote  the  number  of  triangles  (vertices)  about  an  interior  vertex  c.  The  .Vc 
equations  given  by  Eq.  (2.23)  for  the  2A’c  unknowns  {<,+1/2}  and  {u,+i/2}  are  linearly 
dependent.  This  can  be  seen  by  summing  the  equations,  which  produces  the  equation 
for  the  change  in  circulation  about  vertex  c.  The  equation  for  the  change  in  circulation 
at  vertex  c  is  a  linear  combination  of  the  L+i/2's.  which  is  c([ual  to  zero.  Since  we  want 
the  L+i/2's  to  be  linearly  independent,  we  can  set  L+i/2  =  0  for  all  /.  We  still  need 
another  equation  to  determine  the  normal  component  for  the  change  in  velocities  011 
the  triangles. 

Let  us  for  the  moment  write  that  equation  as 


y  C,+  |/2".+  1/2  -  ''• 


(2.211 


Using  Ef|.  (2.23)  with  U  +  1/2  =  0  for  all  /,  we  can  snccessi\'el\'  l■lilllilNlte  e.u  h  fur 

i  —  1 . .Vc  —  1  in  Er|.  (2,24)  until  we  arrive  at  an  ec|uaii(>ii  I'm-  \)  \  ^|  /j.  .^iiice  Hie 


numbering  of  the  triangles  and  vertices  is  arbitrary,  this  expression  is  valid  for  ea<  li 
triangle  i  +  1/2  by  replacing  n  v^.,.1/2  with  n,+i/-.>  and  V  v^^i/j  witli  The  result 

is  that 


^'^i+1/2 


X  (r.- 


+  1 


k{c) 


c*.-+i/2kr-+i  -  t'fci 


•2.4,. 


*■■+1/2 


(vjt+i/2  -  'V.  +  1/2)  • 


/z- 


T-s-i/jIft'+i  -  rrl 


k{c) 


2.4 


r-s- 1 /2 


(2.25) 


Several  alternatives  are  possible  for  Eq.  (2. ‘24).  If  we  conserve  divergence  about  the 
vertex  c,  then 


<^1+1/2  —  l*‘*+i  ’‘‘li' 

6  =  0. 


(2.26) 


The  transformation  R  prescibed  by  Eq.  (2.25)  is  time-reversible,  hence  Eqs.  (2.10)- 
(2.14)  are  also  reversible.  The  entire  algorithm  advances  vertex  positions  and  veloci¬ 
ties  reversibly  while  evolving  the  correct  circulation  about  every  interior  vertex.  This 
technique  is  unique  for  Lagrangian  codes,  which  usually  either  ignore  conservation  of 
circulation  completely  or  conserve  circulation  through  an  iteration  performed  simul¬ 
taneously  with  the  pressure  iteration.  With  this  method  the  circulation  is  conserved 
exactly  regardless  of  whether  the  pressures  have  iterated  to  their  final  values. 


E.  \'iscous  Flows 

Viscosity  modifies  Eq.  (2.9).  so  that  now 


UV  o 

p—  +  Vp  =  fe  +  pV -V. 
at 


(2.27) 


Discretization  of  the  additional  term  in  the  momentum  equation  follows  the  .same  ap¬ 
proach  as  the  discretization  of  the  other  terms.  Since  the  velocity  is  a  triangle-centered 


ijuantity.  we  need  a  discrete  vertex-centered  gradient  operator,  and  a  discrete  triangle- 
centered  divergence  operator.  Employing  the  same  techniques  as  above  we  have 


■dc(V/)<.  =  ^  ^/,+i/2{r.+i  -  r.) 


.4j(V  ■  v)j  =  ^  X  (r.-i-i  -  r.-i)] 


The  Laplacian  is  found  by  taking  the  divergence  of  the  gradient. 

The  finite  difference  equations.  Eqs.  (2.10)  and  (2.14),  can  be  modified  to  account 
for  the  additional  term  in  the  momentum  equation  by 


=  v° - ^(Vp)"  +  —f  +  ^^^(V-v)° 


(2.30) 


^n.k  ^  -  l/2.f  _  A(Vp);'*’-  +  )  ,  (2.31) 

■'  ■'  2pj  •'  2pj  2pj 

These  equations  are  implicit  in  the  velocities,  just  as  the  original  Eqs.  (2.10)-(2.14) 
are.  .\s  in  the  inviscid  case,  we  solve  by  iteration. 

This  algorithm  was  tested  by  calculating  the  spreading  of  a  shear  layer  of  initially 
zero  thickness  given  by 


{(Pr,0),  fori/>yo. 

(0,0),  iory^yo. 

(-Ur,0),  fory<yo- 


(2.32) 


where  y^  is  the  original  location  of  the  vortex  sheet.  The  velocity  distribution  across 


this  laver  evolves  as 


v{x,y,t)  -  .i  Vj.  erf 


(y  -  yo)' 


where  u  =  p/p.  The  width  Ay  of  the  layer  grows  as 


Ay  sj  8(p't)*^"- 


For  the  test  calculation  the  grid  was  initialized  to  center  a  vortex  sheet  in  a  grid 
16  cells  wide  with  an  initial  layer  width  of  zero.  The  two  opposing  streams  had  initially 
constant  velocity  profiles.  The  evolution  of  the  interface  between  the  streams  was 
governed  by  the  same  algorithms  as  the  interior  of  either  fluid,  so  that  no  special 
interface  boundary  condition  was  used.  The  boundary  conditions  on  the  sides  of  the 
computational  region  were  periodic,  and  the  top  and  bottom  had  frce-slip  boundary 
conditions. 

.4t  the  end  of  the  calculation,  the  layer  v/idth  agreed  to  within  numerical  roundoff 
with  the  theory  and  the  layer  extended  over  the  whole  mesh.  The  velocity  profile  for 
each  stream  coincided  with  that  given  by  Eq.  (2.33)  to  within  round-off  error.  The 
[/-components  of  the  velocity  remained  zero,  indicating  that  the  algorithm  was  working 
well  for  the  grid  distortions  presented  by  the  problem. 

F.  Conservation  of  Vertex  Cell  Areas 

Equations  (2.10)-(2.14)  are  implicit  in  the  triangle  velocities  {vj}.  Because  those 
equations  must  be  solved  iteratively  to  produce  a  divergence  free  \elocity  field,  a 
small  residual  error  may  remain.  In  addition,  vertex  velocities  are  derived  frum  the 
divergence-free  triangle  velocities.  In  practice  this  means  that  vertex  cell  areas  may 
not  be  conserved.  Furthermore,  as  the  flow  progresses,  the  triangle  sides  distort.  \et 
at  any  given  time  we  compute  using  straight  triangle  sides,  which  does  not  produce  the 
equivalent  cell  area  about  any  given  vertex.  However,  since  we  know  what  the  triangle 
area  should  be.  it  is  possible  to  at  least  make  a  correction  to  the  known  error.  Our 
approach,  then,  is  to  perform  an  ad  hoc  correction  step  after  all  the  vertices  have  been 
advanced  in  time.  This  correction  step  moves  the  vertices  in  order  to  conserve  vertex 
cell  area.  .After  this  vertex  correction  step,  the  rotator  is  ap[)hc(.l  to  ensure  that  i  lu' 
circulation  has  not  l^een  changed. 


To  expand  or  contract  a  vertex  cell  area,  we  must  expand  or  contract  the  surround¬ 
ing  triangle  areas.  Suppose  we  wish  to  expand  a  triangle  j  with  area  .4^  and  vertex 
coordinates  r,  by  an  amount  SAj.  To  do  this  we  will  move  each  vertex  r,.  an  amount 


j.'>ew 


r.  =  6r,  =  dj  [i  X  (r._, 


(2.3o) 


that  is,  the  vertices  of  the  triangle  are  moved  normally  to  the  opposite  side  a  distance 
prescribed  by  the  triangle  expansion  factor,  dj.  If  dj  is  positive,  the  triangle  area 
increa.ses.  Using  the  vector  definition  for  the  area  of  a  triangle,  we  have 


2SAj  =  2.4"'="'  -  2.4j 

=  [(CT  -  i-r"')  X  (r^  -  r;--)]  •  f  -  [(r.+,  -  r.)  x  (r.  -  r._, )]  •  i 
=  [(r.+i  -  r.)  X  (^r,_i  -  ^r.+i)]  •  c  +  [(<5r,+i  -  ^r,)  x  (r._]  -  r,  +  i )]  •  f  (2. .30) 
+  [(^r,+i  -  ^r.)  X  ((fr,_i  -  <^r.+i)]  •  ~ 

=  s~dj  +  6.4j(/', 

where  s*  is  the  sum  of  the  squares  of  the  sides  of  the  triangle.  This  quadratic  in  the 
expansion  factor,  dj.  can  be  solved  to  yield 


,  -05-  +  +  48.4_,^.4; 

<i,  -  - - -  U-i.l 

The  sign  in  front  of  the  square  toot  was  chosen  to  ensure  dj  has  the  same  sign  as  f.dy. 

We  relate  the  change  in  triangle  area,  to  the  conservation  of  vertex  cell  areas 
through 


<(j)  <w) 

where  the  sum  is  over  the  three  vertices  of  the  triangle.  .4,  is  the  current  area  aliout  ver¬ 
tex  i  and  .4-=  is  the  original  area  about  vertex  I.  Basically,  the  change  in  \ertex  cell  areas 
is  apportioned  to  each  contributing  triangle  according  to  that  triangle's  contrilmt ion 
to  the  vertex  cell  area. 


Although  this  residual  correction  is  a  small  numerical  efTect.  we  have  found  that  it 
improves  the  overall  results  of  a  calculation.  Because  this  algorithm  e.vpands  triangles, 
it  has  potential  for  modelling  other  physical  processes.  In  a  compressible  algorithm 
involving  energy  release  and  Iluid  flows  with  transit  times  which  are  small  compared  to 
the  energy  release  times,  this  algorithm  could  be  used  to  produce  the  reciuired  e.xiKuision 
of  the  vertex  cells. 

G.  Grid  Restructuring 

In  Lagrangian  calculations  the  grid  may  distort  to  the  point  where  grid  restructur¬ 
ing  is  necessary.  The  derivation  of  the  reconnection  and  vertex  addition  and  deletion 
algorithms  are  done  through  the  control  volume  approach  and  the  use  of  triangle  veloc¬ 
ities.  For  all  the  algorithms  used,  the  area-weighted  divergence  and  curl  taken  about 
each  vertex  are  both  identically  conserved  for  grid  reconnections  and  vertex  addition. 

The  accuracy  of  a  general  triangular  mesh  is  diminished  by  large  obtuse  angles 
within  triangles.  With  reconnections,  accuracy  can  be  recovered  by  ensuring  that  large 
obtuse  angles  are  preferentially  eliminated.  There  are  many  ways  of  formulating  a 
reconnection  algorithm.  The  one  we  have  chosen  is  based  on  requirements  for  solving 
the  pressure  Poisson  equation.  The  pressure  Poisson  equation  is  formally  equivalent  to 
that  obtained  by  a  piece-wise  linear  Rayleigh-Ritz-Galerkin  finite  element  procedure 
on  a  triangular  grid.  (See,  for  example  [15].)  Since  we  solve  the  equation  by  iteration, 
we  want  the  iteration  to  converge  as  rapidly  as  possible.  Mathematically,  convergence 
is  assured  if  the  finite  difference  equation  has  a  maximum  principle:  that  is,  all  the 
off-diagonal  terms  are  negative,  the  diagonal  term  is  positive  and  greater  than  or  t'qual 
to  the  absolute  value  of  the  sum  of  the  off-diagonal  terms,  with  strict  inefiuality  for  at 
least  one  equation.  (That  one  equation  typically  involves  boundary  londitions.  Our 
boundary  condition  prescribes  the  integrated  pressure  along  the  upper  Iroundary. ) 


To  see  how  large  angles  affect  the  maximum  principle,  consider  the  difference 
equation  for  vertex  I  of  Fig.  3a.  The  off-diagonal  coefficient  relating  vertex  /  to  vertex 


a  =  --(cot0  -t-cot^"*") 


(2.39) 


where  0+  and  0~  are  the  angles  opposite  the  line  from  the  vertex  j  to  the  vertex  / 
as  shown  in  Fig.  3a.  The  other  off-diagonal  terms  are  determined  in  a  similar  man- 
nerfrom  the  remaining  edges  eminating  from  vertex  1.  The  diagonal  coefficient  is  the 
negative  of  the  sum  of  the  off-diagonal  terms.  For  positive  area  triangles.  0'^  and  0~ 
are  both  between  0  °  and  180  “ .  Hence,  each  term  in  Eq.  (2.39)  is  negative  only  when 
0'^  +  6~  >  180  ° .  since 

sin(0+-f0_) 


2  sin  6+  sin  0~ 


(2.40) 


The  reconnection  algorithm  ensures  that  the  angles  subtended  by  any  given  edge 
sum  to  no  more  than  180  ° .  If  0+  is  greater  than  180  ° ,  the  grid  line  is  reconnected 
as  shown  in  Fig.  3b.  The  new  angles,  and  9'~ ,  must  sum  to  less  than  180  °  since 
{0++0-+0'++0'-  ' )  is  the  sum  of  the  interior  quadrilateral  angles,  which  must  be  3C0  '  . 
By  chosing  the  diagonal  which  divides  the  largest  opposing  angles,  the  reconnection 
algorithm  preferentially  eliminates  large  angles  in  triangles. 

Interface  sides  are  never  allowed  to  reconnect.  In  such  cases  vertex  addition  algo¬ 
rithms  are  needed.  Vertex  addition  algorithms  are  also  needed  where  the  flow  naturally 
depletes  vertices.  For  vertex  addition,  satisfaction  of  conservation  integrals  is  partic¬ 
ularly  simple.  The  vertex  added  at  the  centroid  of  a  triangle  subdivides  that  triangle 
into  three  smaller  triangles.  .4  vertex  added  to  the  midpoint  of  a  side  subdi\ides  the 
two  adjacent  triangles  into  four  smaller  triangles.  If  the  new  triangle  xelocitivs  are  all 
the  same  as  the  velocity  of  the  subdivided  triangles,  all  conservation  laws  arc'  satis¬ 
fied.  Since  the  reconnection  algorithm  is  also  conservative,  siilisequent  reconnect  ion- 
to  other  vertices  ensure  that  the  only  effect  of  the  addition  is  tin  increase  in  rc'solutioii. 


The  case  is  not  as  obvious  for  vertex  deletion.  Reconnections  can  Ije  used  to 
surround  any  interior  vertex  within  a  triangle.  The  vertex  is  then  removed  and  tlie  new 
larger  triangle  given  a  velocity  which  is  the  area-weighted  sum  of  the  old  velocities, 

.4,v,  =  .4,v.  -t-  Aj^j  +  AkWk.  (2.41) 

Such  a  substitution  redistributes  circulation  in  accordance  with  area  coordinates.  Fig¬ 
ure  4  illustrates  the  triangles  before  and  after  vertex  removal.  If  (4  is  the  vorticity 
about  vertex  4  before  it  is  removed,  then  the  vorticity  about  eacli  of  the  other  three 
vertices  is  increased  by  an  amount  ('  given  by 

Cl  =  A,  (4 /A, 

=  AkG/Ai ,  (2.42) 

=  AiQ/Ai 

where 

C(+C2+C3=C4 

since 

A,-  4-  Aj  -t-  Ak  =  A/. 

Therefore,  total  vorticity  is  conserved  and  redistributed  in  a  reasonable  and  natural 
manner. 


III.  Surface  Tension 


I 

I 


I 


.4.  The  Algorithm 

The  surface  tension  at  an  interface  between  two  materials  depends  on  the  curvature 
of  the  interface.  In  the  conventional  numerical  representation  of  surface  tension,  it  is 
cast  into  a  finite-difference  form  by  fitting  vertices  on  the  material  interface  to  some 
parametric  function.  This  function  is  then  used  to  find  an  estimate  of  local  curvature. 
Once  the  curvature  is  known,  a  surface  tension  force  is  evaluated  and  used  to  accelerate 
interface  vertices. 

This  scheme  fails  in  SPLISH  for  two  reasons.  First,  the  interface  vertices  are 
accelerated  directly  by  surface  tension  forces  evaluated  on  the  vertices.  Since  velocities 
are  centered  on  triangles  in  SPLISH,  the  velocity  field  sees  the  effect  of  the  acceleration 
a  half-timestep  later,  unless  a  secondary  calculation  is  made.  .4s  a  result,  the  pressure 
calculated  within  the  droplet  is  inconsistent  with  that  found  from  the  surface  tension 
formula.  Second,  since  the  pressure  gradient  forces  and  surface  tension  forces  are  not 
calculated  in  the  same  manner,  numerical  errors  result  which  grow  with  each  timestep. 

Both  of  these  problems  are  eliminated  by  a  different  formulation  of  surface  tension, 
in  which  a  surface  tension  potential  is  used  to  generate  the  forces.  The  surface  tension 
force  is  formulated  as  a  gradient  of  a  potential  present  only  at  the  surfaces.  With  this 
method,  the  pressure  gradient  forces  are  calculated  in  the  same  manner  and  on  the  same 
grid  as  the  forces  derived  from  the  surface  tension  potential.  Therefore  both  the  surface 
tension  potential  and  the  pressure  are  dynamically  similar,  and  the  physical  pressure 
drop  across  the  interface  must  e.xactly  cancel  the  surface  tension  forces.  Preliminary 
aspects  of  this  work  were  described  by  Fritts  et  al.  [16,  17]. 

The  finite-difference  algorithms  for  surface  tension  are  straightforward.  The  sur¬ 


face  tension  forces  are  included  through  Laplace's  formula  for  the  pressure  jump  across 
an  interface  [18], 


where  p,  is  the  pressure  just  inside  the  droplet  at  tlie  interface,  po  is  tlie  pressure 
just  outside  the  droplet  at  the  interface,  a  is  the  surface  tension  coefficient  associated 
with  the  two  media  which  define  the  interface,  and  R  is  the  radius  of  curvature  in  the 
two-dimensional  plane.  The  radius  of  curvature  is  positive  at  points  on  the  interface 
where  the  droplet  surface  is  convex  (a  circle  is  convex  everywhere)  and  negative  when 
the  surface  is  concave.  These  pressure  jumps  are  included  in  the  Poisson  eejuation  for 
the  pressure.  The  average  pressure,  (p,  +  Po)/2,  is  computed  at  the  interface  vertices. 
From  the  average  pressure  and  the  pressure  jump,  we  can  compute  a  pressure  gradient 
centered  on  triangles,  both  inside  and  outside  the  surface.  This  pressure  gradient  is 
used  in  the  momentum  equation. 

The  radius  of  curvature  is  computed  from  a  parametric  cubic  spline  interpolant  to 
the  interface  vertices.  Past  calculations  of  droplets  oscillating  due  to  surface  tension 
forces  [19.  20]  also  use  cubic  spline  interpolation.  However,  they  divided  the  surface 
into  at  least  four  segments  (the  top,  bottom,  right  and  left  sides  of  the  droplet)  to 
produce  an  interpolant  on  each  segment.  Each  interpolant  was  matched  at  the  joints 
to  produce  an  overall  curve.  The  parametric  interpolant  used  here  does  not  require 
this  special  matching. 

The  parametric  spline  is  produced  in  the  following  manner.  Denote  the  interface 
vertices  by  r,  =  (.r,. !/,),  i  =  l,...jV,  with  r.v  =  rj.  .Also  define  a  pseudo  arc  length 
parameter,  s.  such  that  the  spline  knots  occur  at  the  points 

=  0, 

(;i.2) 

s.  =  s,_i  -I-  |r,  -  r,_i|.  i  =  2 . V, 


We  generate  the  twice  differentiable  periodic  spline  interpolants.  r(  s)  =  (.r( s ).  i/(  ,s) ) 
from  the  data  {-s,},  and  {r.},  i  =  1 . .V,  as  prescribed  by  DeBoor  [21].  The 


curvature  is  then  given  by 


K  =  R-^  = 


22 


where  the  prime  indicates  difTerentiation  with  respect  to  the  parameter  s.  The  sign  of 
R  at  an  interface  vertex,  r,,  is  given  by  the  sign  of  z  ■  [(r.+i  —  r,)  x  (r,_i  —  r,)]. 

We  can  iterate  the  process  if  necessary.  From  the  spline  fit  we  can  generate  new 
values  for  the  {s,}  by  integrating  the  expression  for  arc  length  along  a  parametrically 
prescribed  curve.  For  symmetrically  placed  vertices  on  a  symmetric  droplet,  however, 
we  have  found  the  iteration  on  arc  length  parameter  is  unnecessary. 

The  parametric  spline  fit  is  also  used  for  regridding.  When  the  regridding  algorithm 
calls  for  the  bisection  of  a  triangle  side  which  borders  the  two  media,  a  new  vertex  is 
added  on  the  spline  interpolant  between  the  vertices.  This  is  done  rather  than  bisecting 
the  straight-line  segment,  since  a  straight-line  bisection  introduces  spurious  interface 
oscillations.  Bisecting  the  spline  maintains  a  better  overall  shape  for  the  interface. 

B.  Test  Results 

We  tested  the  algorithm  for  surface  tension  in  SPLISH  using  two  test  problems. 
The  first  test  problem  consists  of  internal  capillary  waves.  In  the  second  test  problem 
we  calculated  the  oscillation  of  a  droplet  due  to  surface  tension.  For  completeness  wo 
also  present  calculations  of  internal  gravity  waves  as  a  test  of  the  overall  hydrodynamic 
algorithms  in  SPLISH. 

1.  Internal  Gravity  and  Capillary  Waves 

The  linear  theory  for  the  small  amplitude  oscillation  of  an  interface  between  two 
fluids,  bounded  above  and  below  by  solid  walls,  gives  the  frequency  as  a  function  of 
wavenumber  A:, 

,2  _  (P  ~  I  j  j 

pcoth  fell  +  p'  coth  kh' 

Here  the  upper  fluid  is  of  depth  h'  and  density  p' .  the  lower  lluid  is  of  de[)tli  li  <in<l 
density  p.  cj  is  the  acceleration  due  to  gravity  and  a  is  the  coetiicient  of  surtace  tension 
for  the  two  media.  Following  the  free-surface  wave  cahmlat  iotis  of  Fritts  and  Boris  dl. 


V 
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we  take  k  =  'IirfX,  A  =  2.0  cm,  h  —  h'  =  1.0  cm.  f>  =  2  g/cc.  and  //  =  1  g/cc.  For  an 
internal  gravity  wave,  we  have  g  =  980  cni'/s  and  <7  =  0  dynes/cm.  For  an  internal 
capillary  wave,  we  have  g  =  0  cm'/sec  and  a  —  30  dyne.s/cm.  These  values  give  a 
period  r  =  2~/^'  —  0.22073  s  for  the  internal  gravity  wave  and  r  =  0.50196  s  for 
the  internal  capillary  wave.  The  amplitude  of  the  oscillation  is  taken  as  A  —  0.0672/i. 
For  this  amplitude  the  free-surface  oscillations  of  Fritts  and  Boris  [1]  showed  negligible 


non-linear  effects.  Figure  5  shows  the  initial  grid  for  the  mesh  size  6s  =  0.125  cm. 


Figure  6  shows  the  wave  period  as  a  function  of  mesh  size  for  the  internal  gravity 
wave  problem.  The  ratio  of  timesteps  for  any  two  calculations  was  the  same  as  that 


for  the  mesh  sizes.  Each  data  point  on  the  curve  is  an  average  over  several  periods  and 
is  accurate  to  three  digits.  If  we  extrapolate  to  zero  mesh  size  using  a  parabolic  least- 
squares  fit.  r  =  To  -h  bSs  +  a{6s)-  to  the  data  points,  we  obtain  tq  =  0.2214.  b  =  0.0726. 
and  a  =  0.1549  for  this  problem.  The  extrapolated  value,  tq,  is  accurate  to  0.3'T.  The 
finite-difference  derivatives  given  in  section  II  are  accurate  to  second  order  in  the  mesli 
size  for  triangular  grids  in  which  the  centroid  of  a  vertex  cell  is  the  vertex  itself.  The 
truncation  error  is  linear  in  the  distance  between  the  vertex  and  the  centroid  of  the 


vertex  cell.  This  truncation  error  can  occur  in  this  problem  for  vertex  cells  near  the 
interface  in  our  discretization  and  hence  the  linear  term  in  6s  in  the  alcove  quadratic 
expression.  This  linear  term  has  a  coefficient  on  the  order  of  the  wave  amplitude  which 
is  the  approximate  distortion  of  the  grid.  The  order  of  convergence  for  the  algorithm 
is  essentially  quadratic  with  a  small  linear  contribution. 

Figure  7  shows  the  wave  period  as  a  function  of  mesh  size  for  the  internal  capillary 
wave  problem.  Here  the  least-squares  fit  to  the  data  gives  7o  =  0.4995.  b  =  0.2198. 
and  a  =  0.0640.  The  extrapolated  period  is  accurate  to  0.0%.  With  surface  tension 
included,  the  convergence  is  primarily  linear  in  the  mesh  size.  The  reduction  in  rate  of 
convergence  is  is  due  to  the  use  of  cubic  splines  to  calculate  the  curxatures.  Flic  cubic 
sidine  curve  itself  is  fourth-order  accurate,  and  theorems  exist  showing  the  seccuid-< u diu' 


accuracv  of  its  second  denvatiVL-s.  Ilou-evcr.  wo  know  of  no  iheoiein  "ivniR  the  accuracv 


of  the  combination  of  derivatives  needed  to  produce  the  curvature  in  Eq.  (d,3j. 


Droplet  Oscillation. 


As  a  futher  test  of  the  algorithm  fur  surface  tension  in  SPLISII.  we  calculated  th 


oscillation  of  a  droplet  due  to  surface  tension.  Rayleigh  [22]  derived  a  linear  theory  for 
small  amplitude  oscillations  on  cylindrical  jets  that  applies  to  the  cylindrical  droplets 
we  are  discussing.  He  concluded  that  when  the  perturbation  is  totally  in  the  plane 
perpendicular  to  the  axis  of  the  cylinder,  the  frequency,  u,'.  for  the  oscillation  is  given 


2  /  3  \ 

=  (n  -  n)  —  . 

per 


where  the  surface  of  the  droplet  is  given  in  polar  coordinates  by 


r  =  a  +  ecos(n^). 


where  p  is  the  density  of  the  jet.  a  is  the  unperturbed  radius  of  the  jet.  and  n  prescribes 
the  mode  of  oscillation  in  the  plane  with  amplitude  c.  For  large  amplitude  oscillations. 
Rayleigh  found  that  the  experimental  frequency  diverged  from  that  predicated  by  the 


linear  theorv.  and  he  attributed  these  differences  to  nonlinear  effects. 


We  have  extended  Rayleigh's  theory  to  include  the  presence  of  an  external  fluid. 
Ecpiation  (3.4)  then  becomes 


■I  =  -  '0 


\p,i  -t-  Pc  ' 


where  pj  is  the  droplet  density  and  p^  is  the  density  of  the  external  fluid. 


The  tests  of  the  surface  tension  algorithm  consisted  of  a  .series  of  calculations  of 


oscillations  initiated  in  the  lowest  oscillating  mode,  u  =  2  in  Eq.  (4. fit.  .Mso.  wv  have 


chosen 


(I  —  0.0 12o  cm 


a  =  30  dynes/cin, 

values  which  are  typical  fur  many  practical  droplet  problems.  W'e  discuss  r<‘sull^  f'ui' 
two  difTorent  sets  of  conditions.  First  we  consider  a  droplet  density  of  2  g/vr  in  a 
background  e.xternal  fluii.l  density  of  1  g/cc.  If  we  use  the  definition  of  the  period  a.' 
2“/w',  Eq.  (3.C)  gives  a  period 

r  =  1.13  X  10"^  s. 

The  second  set  of  conditions  are  for  a  kerosene  droplet,  with  density  0.82  g/cc.  in  a 
background  of  air.  with  density  0.0013  g/cc.  This  second  case,  with  the  650:1  density 
ratio,  is  a  stringent  test  of  the  numerical  approximations. 

Figure  8  is  a  composite  of  frames  from  a  calculation  in  which  e  =  0.2a  =  0.0025  cm 
for  the  2:1  density  ratio  case.  In  this  calculation  there  are  IT  vertices  in  each  direction 
along  the  exterior  boundaries,  12  vertices  on  the  droplet  interface  and  a  total  of  313 
vertices  initially  in  the  calculation.  The  computational  domain  is  0.1  cm  on  a  side.  The 
left  and  right  boundaries  are  periodic  while  the  top  and  bottom  boundaries  are  solid 
walls.  The  timestep  is  bt  =  2.5  x  10“^  s.  The  figures  shows  four  and  a  half  oscillations 
of  the  droplet.  We  can  see  that  as  the  calculation  proceeds,  no  new  vertices  have  been 
added,  but  in  fact  some  have  been  subtracted.  This  was  the  case  because  the  initial 
gridding  was  adequate  to  represent  the  droplet  shape.  From  these  calculations,  the 
period  of  oscillation  is 

rj  i  =  1.35  X  10“^  s. 

Similar  calculations  with  20  vertices  surrounding  the  <lroi)lel  (a  21x21  grid)  shuw  a 


period  of 


r.,0  =  1.33  X  10“  ‘  s. 


for  21  vertices  surrounding  the  droplet  (a  25x25  grid)  we  ha\’e  a  perital  of 


724  =  1-31  X  10“'  S. 


I 


and  for  28  vertices  surrounding  the  droplet  (a  33x33  grid)  the  period  is 

To  8  1 . 2  i  X  10  ^  s . 

In  each  case,  the  period  does  not  change  during  the  calculation.  Figures  0  and  10  'lunv 
the  initial  oscillation  for  the  more  resolved  cases.  For  these  calculations,  it  was  neccs.oary 
to  decrease  the  timestep.  as  discussed  below.  The  time  step  for  the  calculation  with  12 
vertices  surrounding  the  droplet  is  such  that  the  period  cannot  be  resolved  to  better 
than  two  digits.  It  appears  that  the  calculations  are  not  converging  to  the  theoretical 
value,  but  to  a  value  of  1.19  ±  .06s,  based  on  the  graph  of  the  computed  period  as  a 
function  of  mesh  size  shown  in  fig.  11.  The  convergence  is  essentially  linear  as  it  was 
in  the  internal  capillary  wave  test  problem,  but  with  a  numerical  error  of  about  .5.5% 
for  this  calculation. 

Since  the  internal  wave  tests  show  much  better  convergence  properties  for  the  algo¬ 
rithm.  as  do  previous  free-surface  wave  calculations  [1.  2],  than  the  droplet  oscillation 
test  problem,  we  performed  several  other  numerical  tests  on  the  droplet  oscillation 
problem  to  determine  if  the  poorer  convergence  properties  were  due  to  other  numerical 
parameters. 

Firstly,  we  tested  whether  the  presence  of  boundaries  a  finite  distance  away  c.'ould 
alter  the  calculated  period  by  performing  calculations  in  a  larger  domain  of  length 
0.2  cm.  Here  there  were  twice  as  many  vertices  on  the  boundary,  but  still  only  12 
vertices  surrounding  the  droplet  which  was  the  same  size  as  the  droplets  in  the  tests 
described  above.  These  calculations  showed  no  change  in  period,  so  we  conclmle  that 
the  effects  of  periodic  boundaries  and  reflecting  walls  are  negligible. 

It  was  also  important  to  evaluate  the  possible  effects  of  iioiiline.u  it y  m  t  he  s..,:i'  i.ai. 
File  theoretical  value  is  from  a  linear  analysis,  and  the  calculaiiuu  i^  a  full  mwi.  ui''.;; 
calculation.  It  is  possible  that  this  could  account  for  part  of  the  di^crepaia  F- 
this,  we  performed  calculations  with  smaller  amplitudes,  c.  to  'et.'  if  tlaue  ,ii;\ 
difference  in  calculated  period.  The  result  was  that  the  numerical  \alue  of  iIh-  yesH; 


was  the  same  for  e  =  0.01a  =  0.000125  cm  over  the  course  of  two  oscillations  as  it  was 
for  €  =  0.2a.  Our  conclusion  is  that  the  calculations  were  in  a  range  in  which  the  linear 
theory  is  valid. 

We  used  two  diagnostics  to  determine  the  period  of  the  computed  droplet  oscil¬ 
lation.  One  is  the  time  history  of  the  position  of  the  rightmost  vertex  on  the  droplet 
interface,  denoted  by  Xr-  The  other  diagnostic  is  the  quadratic  moment,  defined  by 

<  x"  >=  J  x'dx  dy,  (3.7) 

where  the  integral  is  performed  over  the  triangles  which  define  the  droplet.  Tables  1  - 
IV  give  the  values  of  <  x"  >  and  Xr  as  a  function  of  time  for  the  resolutions  of  12  and 
24  vertices  surrounding  the  droplet.  From  the  maxima  and  minima  in  Tables  I  and  III. 
we  can  determine  the  period  of  oscillation.  It  is  less  well  defined  from  the  values  of  Xr 
in  Tables  II  and  I\h  However,  it  never  differs  by  more  than  two  timesteps  from  that 
gii'en  by  the  moments. 

Finally,  we  examined  the  oscillations  of  a  kerosene  droplet  in  air.  This  calculation 
tests  the  effects  of  the  external  fluid  density  on  the  numerical  convergence  of  tlie  pres¬ 
sure  algorithm  as  well  as  any  role  the  external  fluid  may  have  in  introducing  liigher 
frequency  modes.  Here  the  theoretical  value  of  the  period  is  5.9  x  lO"'*  s.  Using  a  resolu¬ 
tion  of  12  vertices  around  the  droplet  surface,  we  find  a  computed  value  of  7.1  x  I0“  ‘  s. 
The  ratio  between  the  theoretical  and  numerical  results  is  0.S3.  compared  to  a  ratio 
of  0.84  for  the  2;1  density  ratio  calculation  at  the  same  resolution.  Since  changing  the 
density  ratio  from  2;1  to  650:1  did  not  alter  the  relative  error,  we  conclude  that  only 
minor  errors  arise  by  including  of  the  e.xternal  fluid  in  the  calculation. 


C.  Some  Difficulties  and  Limitations 

We  now  believe  the  that  the  inability  of  the  method  to  produce  as  accurate  a 
solution  for  the  droplet  oscillation  test  problem  as  for  the  internal  wave  test  problems  is 
a  combination  of  the  physical  problem  itself  and  the  spline  approximations  to  curvature. 

The  surface  tension  algorithm  discussed  above  suffers  a  basic  problem  in  curve 
fitting.  We  are  trying  to  approximate  an  unknown  continuous  function  by  a  known 
curve  through  a  finite  number  of  points  or  computational  cells.  For  example,  we  are 
trying  to  represent  the  droplet  interface  or  capillary  wave  interface  by  a  spline  fit  to 
a  finite  number  of  points.  Whereas  an  accurate  interpolant  can  be  found  that  goes 
through  a  set  of  points,  it  is  not  always  clear  that  the  other  properties  of  the  curve 
calculated  at  the  points,  e.g.,  the  curvature,  are  well  represented  by  this  interpolant. 
Splines  are  notorious  for  introducing  spurious  oscillations  between  the  points  defining 
them  initially. 

Figure  12  shows  the  curvature  at  each  vertex  around  the  droplet.  The  exact 
curvature  for  the  initial  drop  is  compared  to  the  curvature  produced  by  the  spline 
interpolant  and  to  curvatures  produced  after  one  oscillation  is  completed.  The  initial 
curvature,  defined  by  splines  on  the  interface  vertices,  is  reasonable.  However,  by  the 
end  of  a  cycle,  there  are  spurious  oscillations  even  though  the  curvature  has  the  same 
basic  shape. 

In  the  internal  capillary  wave  problem,  the  range  of  values  for  the  interface  cur\a- 
ture  was  a  factor  of  15  smaller  than  for  the  droplet  oscillation  problem.  .A.s  a  result  the 
interface  curvature  for  the  internal  capillary  wave  is  determined  with  greater  accuracy. 
In  the  droplet  oscillation  problem  where  the  interface  "bends"  more  sharpl}-.  the  si.)liiie 
has  a  greater  difficulty  approximating  the  curvatures  accurately. 

Interpolations  can  cause  other  problems  in  the  calculations.  Our  calculations  have 
shown  that  the  final  result  can  be  affected  by  the  location  of  additional  vertices  used 
to  obtain  a  better  initial  approximation  of  the  droplet  inter  face.  The  grid  initialization 
procedure  involves  two  phases:  a  first  phase  to  generate  a  course  grid,  and  a  'n  wuii 


phase  which  refines  the  grid  produced  in  the  first  phase.  During  the  refinement  phase, 
we  have  two  chices  for  the  location  of  new  interface  vertices.  The  initial  grid  produced 
in  Fig.  8  placed  new  vertices  on  the  droplet  defined  by  the  Rayleigh  oscillation  mode. 
We  could  also  add  the  vertex  on  the  existing  spline  interpolant.  Figure  13  shows  <  .r  ’  > 
as  a  function  of  time  for  the  two  types  of  initialization.  The  curve  labelled  1  is  the 
calculation  in  which  the  additional  vertices  were  placed  on  the  Rayleigh  drop.  The 
curve  labelled  2  placed  the  additional  vertices  on  the  spline  fit.  .-Vfter  one  oscillation, 
the  value  of  <  >  differs  by  8%.  .A.fter  one  oscillation  the  value  of  <  x-  >  on  the 

curve  with  label  1  is  lower  than  the  initial  value  of  <  x*  >  and  the  value  of  <  x’  >  on 
the  curve  with  label  2  is  higher  than  the  initial  value  of  <  x"  >.  Notice  also  that  the 
period,  as  well  as  the  amplitude,  is  affected  by  the  type  of  initialization. 

In  Fig.  8  we  see  that  the  amplitude  of  the  droplet  oscillation  decays  as  a  function 
of  time  even  though  the  period  is  not  changing.  The  damping  rate  is  about  18%  per 
oscillation.  The  shape  of  the  droplet  at  the  end  of  the  calculation  is  notably  different 
than  it  was  at  that  same  place  in  an  earlier  oscillation  cycle.  In  the  ideal  case,  this 
would  not  occur. 

The  decay  of  the  oscillation  is  also  apparent  in  the  moment  <  x'  >  and  the 
variation  in  the  location  of  Xr  from  oscillation  to  oscillation.  It  is  apparent  from  Fig.  1-1 
that  the  <  x^  >  moment  is  dissipating,  and  it  is  apparent  from  Fig.  15  that  the  overall 
shape  of  the  droplet  is  changing.  Energy  associated  with  this  lowest-mode  oscillation 
is  going  into  other  modes,  which  is  reflected  in  the  reduction  of  the  timestep  ro<iuired 
to  keep  the  computations  stable.  In  general,  to  carry  out  these  droplet  oscillation 
calculations  it  was  necessary  to  reduce  the  timestep  to  the  point  where  we  could  ic-wUc 
the  highest  mode  of  oscillation  the  droplet  could  support  at  a  given  resolution.  When 
we  doubled  the  resolution  around  the  droplet,  we  found  that  the  timestep  had  to  Ke 
decreased  by  a  factor  of  about  about  2.8.  This  is  consistent  with  the  analvsis  wliidi 
says  that  since  the  period  is  inversely  proportional  to  where  n  is  the  mode  .<! 


oscillation.  Increasing  the  resolution  of  the  droplet  interface  by  a  factor  of  2  means 
that  the  timestep  must  decrease  by  a  factor  of  2^''-  %  2.8. 


.A.  physical  mechanism  for  the  observed  decay  in  the  n  =  2  normal  mode  oscillation 
is  the  e.xistance  of  a  resonance  between  the  n  =  2  and  the  u  =  3  normal  modes;  that  is 
w'3  =  ±2u.’2.  .a  similar  behavior  in  three  dimensions  has  been  analyzed  by  .N’atarajan  and 
Brown  [23].  In  their  three  dimensional  analysis,  significant  energy  can  be  transferied 
from  one  resonant  mode  to  another  within  ten  oscillations. 

In  summary,  the  total  damping  rate  for  the  droplet  oscillation  calculation  is  roughly 
18%  for  the  n  =  2  mode.  Much  of  the  energy  from  this  mode  is  transferred  to  higher 
harmonics,  as  evidenced  by  the  calculated  droplet  shapes,  curvatures  and  the  numerical 
timestep  limitations.  The  difference  in  initialization  procedures  alone  produced  an  8% 
change  in  amplitude.  Since  the  total  numerical  error  is  5. .5%,  we  conclude  that  the 
majority  of  this  error  arises  from  the  inability  of  the  spline  fit  to  appro.ximate  large 
curvatures  accurately.  This  error  is  large  enough  to  mask  other  error  contributions, 
so  that  we  cannot  evaluate  additional  error  terms  other  than  to  indicate  that  they  are 


apparently  much  smaller  than  that  due  to  the  spline  fit. 

However,  despite  all  the  problems  with  spline  fits,  we  found  that  they  provided  a 
good  way  to  calculate  curvature.  In  the  search  for  better  curvatures,  we  have  also 


tried  other  methods,  none  of  which  produced  better  results.  We  enumerate  these 
attempts  both  for  completeness  and  to  emphasize  that  better  numerical  approximations 
for  curvature  are  still  needed  to  permit  more  accurate  calculations  of  surface  tension. 

(1)  We  averaged  the  curvature  between  the  vertices,  such  that 


^  i  +  I  ^1  —  1 


K{s)cU. 


The  results  were  found  to  depend  sensitively  how  how  rhe  integial  wa^  ,ii  t>i,dl\ 
performed.  However,  integration  produced  results  whiidi  were  no  betti'r  than  tin' 
[)ointwise  curvatures  discussed  above. 


(2)  We  smoothed  the  curvatures  K,  with  a  least  squares  linear  sihine.  I  his  nieiliud 
worked  well  for  one  period,  but  the  method  failed  on  subsefiueni  ijsrillatiun-. 

(3)  We  used  a  circular  arc  to  calculate  the  curvatures,  circle  was  placed  tiirwuuli 
the  three  adjoining  vertices.  The  radius  of  that  circle  was  used  as  the  radius 
of  curvature  for  the  interface  at  the  center  verte.x  of  the  three  vertices.  I  his 
method  did  not  work  at  all.  The  droplet  interface  distorted  wildly  w'ithin  the  tirst 
oscillation. 

(4)  We  used  splines  under  tension.  This  approach  introduced  a  free  parameter  winch 
could  not  be  consistently  determined. 

(5)  Based  on  the  experiences  of  Foote  [20],  we  tried  producing  an  interpolant  through 
every  other  vertex  and  averaging  the  result.  The  motivation  was  that  fewer  points 
could  introduce  fewer  oscillations,  and  that  averaging  the  interpolants  could  damp 
the  oscillations.  This  produced  poor  results.  The  calculation  is  really  the  average 
of  two  calculations  with  half  the  original  accuracy. 

(6)  We  considered  but  did  not  implement  nonlinear  splines  [24].  .■\lthough  these  splines 
produce  differentiable  curvatures,  there  is  no  guarantee  that  there  exists  such  a 
spline  through  a  given  set  of  points  and.  if  such  a  spline  does  exist,  there  is  no 
guarantee  that  it  is  unique. 

(7)  We  considered  several  methods  for  calculating  an  interpolant  based  on  the  Rayleigli 
modes.  The  high  mode  oscillations  could  then  be  eliminated.  .\’one  of  the 
schemes  we  considered  gave  better  results  than  the  spline  interpolant.  and  all 
introduced  arbitrary  parameters  into  the  calculation.  These  parameters  ewuli!  in- 
well-determined  for  a  particular  known  shape,  but  eould  lea  In'  fii'termn:'''i  1"',  a 
general  unknown  shape. 


IV.  Incompressible  Flow  about  a  Droplet 

In  this  section  we  present  some  preliminary  calculations  of  forced,  asymmetric 
drop  oscillations  induced  by  flow  around  a  droplet.  These  calculations  include  both 
the  effects  of  viscosity  and  surface  tension.  The  capability  of  studying  such  flows  for 
highly  viscous  droplets  in  shear  flows,  in  two  and  eventually  in  three  dimensions,  is  the 
motivation  for  developing  the  viscosity  and  surface  tension  algorithms. 

The  initial  conditions  we  used  specified  an  initially  steady-state  potential  flow 
about  a  periodic  series  of  cylinders.  .A.gain,  the  boundary  conditions  on  the  left  and 
right  sides  are  periodic,  and  the  upper  and  lower  boundary  conditions  are  reflecting 
walls.  Initially,  a  perfectly  circular  droplet  is  at  rest  in  a  background  flow.  A  physical 
situation  modelled  by  such  an  initialization  might  occur  if  the  flow  velocity  were  ramped 
up  to  its  final  value  before  any  significant  structure  could  develop  in  the  flow,  and  before 
the  droplet  could  pick  up  any  substantial  velocity.  Basically,  it  is  a  smooth  start  for  the 
calculation.  Previously  we  had  performed  calculations  which  began  with  an  impulsive 
start,  but  found  that  as  a  result  there  was  a  large  amount  of  momentum  transferred 
across  the  droplet  interface  early  in  the  calculation. 

The  calculations  presented  here  model  the  forced  fluid  flow  due  to  a  fast  air  stream 
about  an  initially  stationary  kerosene  droplet.  The  physical  parameters,  given  in  Ta¬ 
ble  V.  are  appropriate  for  a  combustor  environment.  .A  total  of  .309  vertices  were  used 
to  initialize  the  problem,  with  12  vertices  at  the  droplet  interface.  Figure  IG  fulluws 
the  evolution  of  pathlines  in  the  internal  and  e.Kternal  flow  fields  through  a  >ei;e>  ot 
timesteps.  For  an  air  velocity  of  100  m/s  and  a  droplet  radius  of  125  microns,  ilie  l  ui- 
responding  Reynolds  number  is  roughly  1600.  The  pathlines  are  detined  by  tin'  [i.iiiis 
of  vertices  over  five  timesteps.  By  the  last  frame  of  Fig.  IG.  the  fluid  originally  tw  -  In' 
left  of  the  droplet  has  progressed  through  the  mesh  and  interacted  with  tli<.'  I'acf  .e  :  in' 
( ne.Kt )  droplet. 

The  first  clear  indication  of  the  development  of  the  recirculation  reiriun  -'•'•n 
in  the  fourth  frame  of  Fig.  16,  which  shows  a  pair  of  counter-rotating  vortua  -.  hi.' 


recirculation  zone  continues  to  develop  throughout  the  calculation,  although  at  limt-s 
the  vortex  pair  is  not  as  evid'.'Ut  due  to  the  deletion  and  addition  of  vertices,  wiiii  li 
interrupts  the  continuity  of  the  pathlines.  By  the  last  frame,  another  pair  of  vortic  es 
is  forming  near  the  droplet,  and  the  original  pair  has  been  shed.  The  leading  f.ic  e  uf 
the  droplet  is  now  quite  distorted,  and  the  droplet  is  about  to  enter  the  wake  ,,f  the 
preceding  droplet. 

Distortions  in  the  face  of  the  droplet  are  evident  in  at  least  the  seventh  frame. 
These  distortions  occur  because  the  curvature  has  increased  and  the  streamlines  of 
the  streamlines  in  the  external  flow  are  condensed  by  the  approaching  wake.  The 
internal  velocities  are  small  compared  to  the  external  flow  rates  and  therefore  cannot 
be  distinguished  as  pathlines.  However,  indication  of  the  (small)  internal  recirculation 
may  be  obtained  by  comparing  internal  vertex  positions  at  various  timeteps. 

Figure  17  shows  the  grid  at  times  in  the  calculation  corresponding  to  those  in 
Fig  16.  During  the  course  of  the  calculation,  a  great  deal  of  vertex  addition  and  delei  ion 
has  occured.  Vertex  addition,  however,  is  most  noticeable  in  the  wake  of  the  droplet 
and  around  the  droplet  interface.  Whereas  there  were  300  vertices  at  the  beginning  of 
the  calculation,  there  are  450  at  the  end. 

Figure  18  shows  the  pathlines  for  a  simulation  with  the  air  speed  increased  'o  120 
m/s.  corresponding  to  a  Reynolds  number  of  2000.  The  fluid  now  completely  i)asses 
through  the  mesh.  The  fluid  initially  near  the  droplet  has  completely  passed  the  next 
droplet  by  the  time  the  calculation  was  terminated.  The  initial  flow  about  the  droplet 
is  similar  to  that  shown  above,  except  for  a  more  pronounce<l  flattening  at  ilie  fa>:e 
of  the  droplet  due  to  the  higher  flow  speed.  The  wake  dev(doi:i>  in  min:li  il.o  '.nuc 
manner,  but  it  now  interacts  strongly  with  the  flow  at  the  forward  .'t agnation  : 
the  droplet.  Oscillations  in  the  flow  due  to  the  wake  are  traii-'iniited  to  tlie  ;'a\'.,ii'i 
face  of  the  droplet  and  give  rise  to  fairly  large  perttirbatioii'. 

.\s  seen  in  Fig.  19.  the  computational  grid  needs  further  relinement  at  tl:>  'nne 
because  the  perturbations  cannot  be  resolved  by  the  limits  set  on  ininimnin  'li.ii.gi'' 


size  originally  chosen  for  the  calculation.  A  sign  that  the  calculation  is  un<;ler-ii  >ul\ isi 
ij  that  one  of  the  crests  of  the  surface  wave  is  spanned  by  a  single  triangle,  a  situation 
which  allows  no  communication  of  that  surface  fluid  with  the  interior  of  the  droplet.  In 
order  to  continue  the  simulation,  better  resolution  must  be  obtained  about  the  droplet 
surface.  .Another  algorithm  is  currently  being  included  to  allow  higher  resolution  near 
points  of  large  curvature  at  material  interfaces. 

V.  Summary  and  Conclusions 

This  paper  presented  the  current  algorithms  included  in  the  code  SPLISH.  a  two- 
dimensional  Cartesian  Lagrangian  treatm.ent  of  incompressible  flows  with  a  dynamically 
restructuring  grid.  The  new  rotator  algorithm  is  an  improvement  on  the  one  previously 
used  for  conserving  circulation.  The  residual  algorithm  ensures  conservation  of  the  area 
of  cells.  These  algorithms  together  with  the  original  SPLISH  framework  constitute  an 
e.xtremely  flexible  code  for  calculating  incompressible  flows  in  highly  distorted  geome- 
tries  or  with  obstacles  in  the  flow. 

New  algorithms  for  modelling  the  physical  effects  of  viscosity  and  surface  tension 
have  been  added.  Whereas  adding  the  viscosity  algorithm  was  relatively  straightfor¬ 
ward.  adding  surface  tension  caused  a  number  of  numerical  problems.  Detailed  iMUe.ii- 
marks  of  the  final  algorithm  selected  were  presented  using  internal  capillary  wa\i  ami 
a  Rayleigh  oscillating  droplet  as  test  problems.  The  surface  tension  algorithm.  has«‘d 
on  spline  fits  to  determine  curvature,  allowed  the  droplet  to  oscillate  many  tinm-  and, 
still  maintain  a  constant  period.  The  numerical  tests  on  the  iiit'u  nal  capiilar\  .,\i  ' 
indicate  that  the  surface  tension  algorithm  produced  a  conv.u'U'nce  rate  wi.r  i,  ha- 
ear  in  the  mesh  size,  whereas  the  basic  hydrodynamic  algoritimm  ao'  iiuadr.iti'  a.  ' 
mesh  size  for  ideal  meshes.  The  droplet  oscillation  te.^t  proihmii  ,a.,,.,.'.  , 

'oine  difficult ies  with  the  spline  fits  for  curvature  when  the  uitctface  iiei.im.  -  _ 

distorted. 
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Previous  numerical  calculations  of  oscillating  spherical  droplets  with  sm  hwa  ten¬ 
sion  and  viscosity  using  a  marker-and-cell  method  showed  only  one  oscillati.ni  of  a 
water  droplet  in  air  [20].  and  thus  did  not  give  any  information  about  the  subse<iueut 
behavior  of  the  mode  amplitudes.  These  calculations  used  2.5  times  the  resolution  of 
our  most  resolved  calculations.  Their  calculated  period  differed  from  the  theoretical 
period  by  69^  compared  to  our  12%  for  a  similar  initial  deformation.  Their  \  iscou.s 
calculations  failed  to  damp  as  quickly  as  required  by  theory  which  may  indicate  that 
they  suffer  from  a  similar  problem  of  approximating  curvatures. 

We  presented  calculations  showing  how  a  kerosene  droplet  deforms  and  sheds  vor¬ 
tices  in  the  wake  of  a  shear  flow.  Calculations  of  fluid  flow  in  and  around  fuel  droplets 
are  important  in  the  study  of  spray  combustors.  The  flow  patterns  influence  droplet 
breakup,  evaporation  and  burning  rates. 
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Table  II 


0.0000 

0.62.50E-03 


0.TIT5E-02 

0.2G00E-02 

0..327oE-02 

0.3900E-02 

0.4550E-02 


0.5825E-02 

0.6525E-02 

0.7150E-02 

0.7775E-02 

0,34o0E-02 

0.1J0.70E-02 

U.'„)72.3E-02 

U.1u:37E-01 

ii.l  lOOE-01 

U.1167E-01 


10  X  16 


last  period 


.00127.5 


.0013 


.0013 


.00125 


.001325 


.0013 


0.I500E-01 

0.9974E-02 


0.1325E-02  .001325  0.1486E-01 


0.1046E-01 

0.1463E-01 

0.1064E-01 

0.1423E-01 

0.1078E-01 


0.5225E-02  .001325  0.1400E-01 


0.1104E-01 

0.1392E-01 

0.1131E-01 

0.1378E-01 

0.1140E-01 


.001275  0.1355E-01 


0.1146E-01 

0.1335E-01 

U.1159E-01 

0.1332E-01 


.00131 


0.0000 

0.6.594E-03 

0.1306E-02 

0.1966E-02 


0.3.513E-07 
0.1628E-0T 
0. 3420  E- 07 
0.1637E-07 


Table  IV 
24  X  24  grid 


time 

last  period 

Xr 

0.0000 

0.1500E-01 

0.6594E-03 

0.1038E-01 

1I.12.50E-02 

.0012.5 

0.1503E-01 

().lIJ3lE-02 

0.1015E-01 
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Table  V 


density  of  kerosene 

O.S'2  g/cc 

density  of  air 

0.0013  g/cc 

surface  tension  (STP) 

•30  dynes/cm 

viscosity  of  kerosene 

1.8  centipoise 

\  i3<:osity  of  air 

0.018  centipoise 

<iir  velocity 

100  or  120  m/s 

initial  droplet  velocity 

0.0  m/s 

droplet  radius 

12.5  microns 

Figure  Captions 


Figure  1.  A  section  of  a  triangular  grid  showing  a)  a  material  interface,  b'  a  vertex- 
cell. 

Figure  2.  .A  test  problem  for  conservation  of  circulation,  a)  The  initial  flow  pattern, 
b)  The  velocities  after  a  half-time  step,  c)  The  velocities  after  the  old  rot.itur 
operator  is  applied. 

Figure  3.  a)  Definition  of  the  angles  0'^  and  0~  for  the  diagonal  line  drawn  from  j  to  /. 
bj  The  angles  and  d'~  formed  by  connecting  the  other  quadriliateral  di.i'.:im,ii. 

Figure  4.  a)  \’ertex  4  isolated  within  a  larger  triangle  before  its  removal.  Ij)  The 
larger  triangle  remaining  after  deletion  of  vertex  4  and  three  associated  sides  and 
triangles. 

Figure  5.  The  initial  grid  for  the  internal  wave  test  problems. 

Figure  6.  The  period  r  as  a  function  of  mesh  size  for  the  internal  gravity  wave  test 
problem. 

Figure  7.  The  period  r  as  a  function  of  mesh  size  for  the  internal  capillary  wave  rest 
problem. 

Figure  8.  .A  composite  of  frames  from  a  calculation  of  an  n  —  2  normal  moi.a,' 
droplet  oscillation  with  12  vertices  around  the  droplet:  =  1  g/'ce.  p,/  =  a,  cc. 

'T  =  30  dynes/cm.  a  —  0.0125  cm.  Each  frame  is  u.l  x  0.1  cm". 

Figure  9.  .A  composite  of  frames  from  a  calculation  of  an  n  =  2  normal  mode  dre'ide! 
oscillation  with  24  vertices  around  the  droplet.  Same  conditions  as  in  Fig.  7. 

Figure  10.  .A  composite  of  frames  from  a  calculation  of  an  n  =  2  normal  mode  d.iuplet 
oscillation  with  28  vertices  around  the  droplet.  Same  conditions  as  in  Fig.  7. 

Figure  11.  The  period  as  a  function  of  mesh  size  for  the  droplet  oscillation  proMvni. 

Figure  12.  Curvature  as  a  function  of  vertex  index  around  the  chop  in  the  lo  '  iii 
ivdculation.  tl)  Exact  solution:  (2)  initial  spline:  i3i  after  icic  wscillaiiuii. 


Figure  13.  The  moment  <  ,r-  >  as  a  function  of  time  for  two  initializatiuns:  '  1  .ui 

initial  vertices  on  the  Ra\ieigh  drop:  f2)  initial  reliniiig  \t;rt  icvr,  (ui  ihi- 

iiiterpolant. 

Figure  14.  Tlie  moment  <  x-  >  as  a  function  of  time  in  the  16  x  hi  calculation. 

Figure  1.5.  The  position  of  the  rightmost  vertex.  Xp.  as  a  function  of  time  in  the  l(i  x 

16  calculation. 

Figure  16.  Pathlines  from  a  calculation  of  air  flowing  past  a  deforming,  viscous  kerosene 
droplet.  Surface  tension  forces  are  included  at  the  material  interface.  Heads  id' 
pathlines  are  the  current  vertex  positions  and  the  tails  are  made  up  of  the  previou^ 
five  positions.  The  flow  speed  is  100  m/s  and  Re  s;  1600. 

Figure  17.  Frames  showing  the  triangular  grid  at  the  same  times  as  shown  for  the 
pathlines  in  Fig.  14. 

Figure  18.  Pathlines  from  a  calculation  similar  to  that  shown  in  Fig.  14.  but  with  a 
flow  velocity  of  120  m/s  and  Re  a:  2000. 

Figure  19.  Frames  showing  the  triangular  grid  at  the  same  times  as  shown  for  the 
pathlines  in  Fig.  16. 
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Consider  an  inviscid  shear  flow  on  the  grid  shown  in  Fig.  2a.  Triangles  above 
y  =  0  have  a  velocity  =  -1.  and  those  below  have  a  velocity  Cj.  =  +1.  If  after 
one  step  the  vertices  have  moved  as  in  Fig.  2b,  conservation  of  circulation  through 
Ec[.  (2.18)  imparts  a  y-component  to  the  velocities  for  those  triangles  bordering  the 
shear.  .Although  the  circulation  integral  about  each  verte.K  in  the  grid  is  conserved,  the 
flow  is  now  no  longer  independent  of  y. 

To  obtain  a  better  formulation  of  the  transformation  R  we  must  consider  Eq.  (2.17) 
more  carefully.  Since  Eq.  (2.17)  is  linear  in  the  unknowns  {v,}.  we  can  obtain  the  change 
in  triangle  velocities  by  considering  the  change  produced  by  the  movement  of  a  single 
verte.x  c.  with  coordinates  Fc,  and  sum  the  resultant  expression  over  all  vertices.  It  is 
reasonable  to  assume  that  the  rotator  should  change  only  the  velocities  of  the  triangles 
which  have  c  as  a  verte.x.  .4s  a  result,  conservation  of  circulation  gives 

v,+  i/2'(r"-r,+  i  )  +  v,_i/2-(r.-i  -r")  =  v, 4.1/2 -(r^- r.+i ) +  v,_i/2  •  (r,_i  -  r°)  (2.19) 

for  each  vertex  i  about  c.  We  have  used  r;  =  r'*  =  r°  for  those  vertices  which  are 
stationary.  If  only  vertex  c  moves,  the  cell  area  at  vertex  c  Is  constant,  so  that  vortlclty 
is  conserved  about  vertex  c  as  well.  How’ever.  at  all  neighboring  vertices,  circulation, 
not  vorticity.  is  conserved.  By  introducing  the  notation 

4.1/2  =  V. 4.1/2  -  V, 4.1/2 


and 


Eq.  (2.19)  may  be  rewritten  as 


Svr  =  r" 


4.1/2  •  (r"  -  r.+i)  +  (3v,_i/2  •  (r,_i  -  i-;*)  =  (v,_i/2  -  v, 4.1/2)  ■  6i\. 
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